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INTRODUCTION: 


A  tracked  free-hand  ultrasound  system  is  ideal  for  guiding  many  radiotherapy  procedures, 
as  it  can  be  performed  in  the  treatment  room.  In  partial  breast  irradiation,  the  lumpectomy 
cavity  should  be  localized  during  the  treatment  course.  While  many  structures  in  the 
breast  look  similar  to  the  lumpectomy  cavity  in  ultrasonography,  the  scar  tissue  around 
the  cavity  is  hard  and  can  be  visualized  by  ultrasound  elastography.  To  be  clinically 
successful,  the  elastography  method  needs  to  be  robust  to  the  sources  of  decorrelation 
between  ultrasound  images,  specifically  fluid  motions  inside  the  cavity,  change  of  the 
appearance  of  speckles  caused  by  compression  and  physiologic  motions,  and  out-of-plane 
motion  of  the  probe.  In  previous  annual  report,  we  presented  a  novel  elastography 
technique  that  was  based  on  analytic  minimization  of  a  regularized  cost  function.  The 
cost  function  incorporates  similarity  of  RF  data  intensity  and  displacement  continuity, 
making  the  method  robust  to  decorrelation  noise  present  throughout  the  image.  We 
exploited  techniques  from  robust  statistics  to  make  the  method  resistant  to  large 
decorrelations  caused  by  sources  such  as  fluid  motion.  The  analytic  displacement 
estimation  worked  in  real-time.  The  tracked  data,  used  for  targeting  the  irradiation,  was 
also  exploited  for  discarding  frames  with  excessive  out-of-plane  motion.  Since  then,  we 
have  extended  the  ID  method  to  2D,  meaning  that  we  can  calculate  tissue  displacements 
more  completely.  The  2D  method  follows  the  path  of  the  ID  method  in  optimizing 
regularized  cost  functions  and  is  also  real-time.  We  also  introduced  Kalman  filtering  for 
calculating  strain  images.  The  Kalman  filter  allows  estimating  low  variance  strain 
images,  while  it  doesn’t  eliminate  boundaries  by  over- smoothing  the  strain  images.  A 
major  contribution  has  also  been  made  in  utilizing  multiple  ultrasound  images  to 
calculate  the  elasticity  image.  Displacement  estimation  is  an  essential  step  for  ultrasound 
elastography  and  numerous  techniques  have  been  proposed  to  improve  its  quality  using 
two  frames  of  ultrasound  RF  data.  This  paper  introduces  a  technique  for  calculating  a 
displacement  field  from  three  (or  multiple)  frames  of  ultrasound  RF  data.  To  calculate  a 
displacement  field  using  three  images,  we  first  derive  constraints  on  variations  of  the 
displacement  field  with  time  using  mechanics  of  materials.  These  constraints  are  then 
used  to  generate  a  regularized  cost  function  that  incorporates  amplitude  similarity  of  three 
ultrasound  images  and  displacement  continuity.  We  optimize  the  cost  function  in  an 
expectation  maximization  (EM)  framework.  Iteratively  reweighted  least  squares  (IRLS) 
is  used  to  minimize  the  effect  of  outliers.  We  show  that,  compared  to  using  two  images, 
the  new  algorithm  reduces  the  noise  and  eliminates  ambiguities  in  displacement 
estimation.  The  displacement  field  is  used  to  generate  strain  images  for  quasi-static 
elastography. 


BODY: 

We  are  very  excited  to  report  the  novel  achievements  of  this  research  effort.  The  detailed 
Statement  of  Work  tasks  and  a  description  below  each  task  is  provided  next.  The  SOW 
tasks  are  in  blue.  A  description  of  the  research  effort  follows  each  SOW  item  in  black. 
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1.  Obtain  ultrasound  (US),  US  elastography  (USE)  and  CT  scans  before  the 
start  of  the  radiotherapy. 

a.  Develop  a  tracked  US  system  for  data  collection  (month  1).  We  will 
use  Polaris  optical  tracker,  already  available  in  the  ERC-CISST  for 
tracking.  Polaris  markers  will  be  attached  to  the  US  probe  and  the  US 
probe  will  be  calibrated  with  respect  to  the  trackers. 

Polaris  optical  tracker  gives  very  accurate  displacement  measurements. 
However,  they  require  line  of  sight  (as  they  are  optical)  and  also  they  take 
longer  time  to  set  up  in  a  CT  imaging  room.  Because  of  the  time  constraints 
(we  had  a  short  amount  of  time  to  set  up  the  tracking  device  in  the  CT  room 
before  the  patient  shows  up),  we  instead  used  magnetic  trackers.  Magnetic 
trackers  are  not  as  accurate  as  Polaris  and  suffer  from  distortion  of  magnetic 
field  due  to  presence  of  metals.  However,  they  do  not  require  line  of  sight  and 
also  easier  to  set  up  and  also  are  inexpensive. 


b.  Collect  US  data  from  patient  before  the  PBI  treatment  at  the  same 
time  that  CT  is  collected  (months  2-14). 

We  have  collected  data  from  more  than  15  patients  so  far.  All  data  was 
acquired  in  the  CT  room  while  the  patient  was  lying  down  on  the  CT  bed.  The 
patient  did  not  move  after  ultrasound  data  was  acquired  until  her  CT  scan  was 
finished.  We  have  acquired  both  3D  ultrasound  data  (freehand  data,  no  3D 
probe  is  used)  and  palpation  data  (for  elastography).  More  details  of  the 
system  are  in  the  attached  paper  (MICCAI  conference,  34%  acceptance  rate, 
MedLine  and  ISI  indexed  and  listed).  I  received  travel  award  from  MICCAI 
(only  7  students  received  this  award)  to  present  this  work. 

Our  approach  introduces  minimal  divergence  from  the  original  workflow  of 
PBI  treatment.  We  have  an  approved  institutional  review  board  (IRB)  protocol 
to  obtain  B-mode  and  strain  images  from  patients  who  undergo  lumpectomy. 
We  acquire  the  ultrasound  data  when  the  patients  return  for  the  CT  scan  four 
weeks  after  the  surgery.  The  data  includes  tracked  B-mode  images  scanned 
over  the  lumpectomy  bed,  tracked  real-time  strain  images,  and  RF  dara 
synchronized  with  the  tracking  information  for  off-line  processing.  We  have 
devised  a  data  collection  system  for  this  purpose  shown  in  Figure  1 . 

Before  acquiring  the  CT  scan,  four  CT-compatible  fiducials  are  placed  around 
the  scar  of  the  surgery  on  patient’s  breast.  These  fiducials  are  additional  to  the 
ones  that  are  commonly  placed  on  patient’s  chest  and  the  resting  foam  bed. 
The  foam  bed,  shown  in  Figure  1  maintains  the  configuration  of  patient’s 
body  during  the  treatment.  The  extra  fiducials  are  used  to  locally  align  the  CT 
and  ultrasound  data,  and  the  regular  ones  are  used  for  targeting  at  the  time  of 
irradiation. 
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Figure  1.  The  Polaris  camera  (NDI,  Waterloo,  Canada)  tracks  the  3D 
orientation  of  the  ultrasound  probe  by  detecting  the  4  bright  spheres  attached 
to  the  ultrasound  transducer. 


The  patient  remains  in  the  same  resting  position  after  the  CT  scan.  As  the  first 
step,  the  sonographer  digitizes  the  location  of  the  fiducials  with  a  calibrated 
stylus  by  simply  touching  each  of  the  four  fiducials  with  the  stylus  in  a  certain 
order  prescribed  in  our  protocol.  Afterwards,  the  sonographer  sweeps  the 
ultrasound  transducer  over  the  lumpectomy  bed  collecting  at  least  500  B- 
mode  images.  Then  within  a  few  seconds,  the  program  constructs  the 
ultrasound  volume  from  the  collected  images.  Next,  the  sonograhper  obtains 
individual  tracked  strain  images  from  the  areas  of  interest  by  gently  moving 
the  ultrasound  transducer  up  and  down.  The  visualization  program 
demonstrates  the  flying  strain  image  over  the  ultrasound  volume  (Figure  2). 
Depending  on  each  case,  we  acquire  a  minimum  of  200  strain  images.  Lastly, 
we  record  sequences  of  RF  signals  as  the  sonographer  compresses  and 
decompresses  the  tissue.  The  RF  data  synchronized  with  the  tracking 
information  facilitates  more  elaborate  strain  imaging  techniques. 
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Figure  2.  The  visualization  software  shows  3D  ultrasound  and  strain  images. 

c.  Register  US  to  the  CT  (months  2-14). 

We  advanced  significantly  and  proposed  exciting  new  avenues  for  further 
work  in  developing  novel  techniques  for  elastography,  which  is  the  major 
contribution  of  this  research  proposal.  Therefore,  we  have  not  had  enough 
time  to  focus  on  this  task.  However,  we  have  performed  registration  of  US 
and  CT  using  the  tracking  information  (obtained  from  the  Polaris  camera). 
The  registration  is  included  in  the  data  acquisition  software,  which  enabled 
fast  and  smooth  data  collection  and  visualization.  The  application  interfaces 
with  the  ultrasound  machine  and  the  tracker  (optical  or  magnetic), 
synchronizes,  records,  and  visualizes  the  data.  It  is  implemented  in  C++  with  a 
multi-threaded  scheme  for  better  performance.  We  have  also  developed  a 
simple  and  intuitive  graphical  user  interface  (GUI). 

B-mode  images  and  transducer  locations  are  continuously  stored  in  two 
circular  buffers  in  separate  threads  with  the  highest  rate  possible.  A  time- 
stamp  with  microsecond  accuracy  is  also  recorded  along  with  the  captured 
data.  A  constant  delay  calculated  off-line  is  applied  to  synchronize  the  stream 
of  data.  Acquisition  of  B-mode  images  is  fast,  and  the  region  of  interest  can  be 
scanned  within  a  few  seconds.  The  maximum  rate  of  storing  data  is  bounded 
with  the  minimum  of  frame  rate  of  the  B-mode  images  and  the  tracker. 
However,  the  user  can  set  the  rate  of  data  storage  to  a  lower  rate  through  the 
GUI  if  needed.  Due  to  the  high  data  rate  and  large  frame  sizes  captured  per 
second  (typically  10  to  30  frames  per  second),  immediate  saving  to  the  hard 
drive  may  result  in  unexpected  delays  or  loss  of  frames.  Therefore,  the  data  is 
first  stored  in  the  random  access  memory  (RAM)  of  the  computer  and  stored 
in  hard  disk  only  when  the  data  collection  is  stopped. 
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Once  the  B-mode  scan  is  acquired,  the  program  constructs  a  3D-US  volume. 
The  volume  is  presented  as  three  orthogonal  slices  that  can  be  freely  translated 
and  rotated  in  space  as  shown  in  Figure  2.  VTK  (Kiteware  Inc.)  is  mainly  used 
for  the  visualization  and  reslicing  tasks  (Figure  3).  The  real-time  strain  images 
are  captured  in  the  same  manner  except  that  the  strain  plane  is  overlaid  on  3D- 
US.  In  addition,  the  program  can  sends  commands  to  ultrasound  machine  for 
capturing  sequences  of  RF  data.  When  the  ultrasound  machine  receives  the 
command,  captures  and  saves  one  sequence  of  RF  frames,  while  our  program 
records  the  tracking  information. 


(a)  Resliced  CT  (b)  Magnified  CT  (c)  B-mode  image 


Figure  3.  Registration  of  CT  and  ultrasound,  (a)  Shows  a  CT  image,  sliced  in 
3D  (using  interpolation)  to  match  the  ultrasound  image,  whose  orientation  is 
known  from  the  tracking  information.  In  (b),  the  square  region  of  (a)  is 
magnified,  (c)  shows  the  corresponding  ultrasound  image. 

In  addition,  we  performed  some  efforts  to  improve  the  tracking-based 
registration  results  of  Figure  3  by  generating  a  pseudo-ultrasound  image  from 
the  CT  volume  and  registering  the  CT  and  ultrasound  through  the  pseudo¬ 
ultrasound,  but  the  preliminary  results  were  not  satisfactory,  and  we  decided 
to  use  the  tracking-based  registration  results  of  Figure  3. 


d.  Compare  different  combination  of  US,  USE  and  CT  for  delineation 
(months  3-18). 

We  have  done  preliminary  studies  on  comparing  the  performance  of  the 
operator  with  different  imaging  modalities.  While  it  is  intuitive  that  adding 
ultrasound  and  elastography  enhance  the  delineation  of  the  lumpectomy  bed, 
we  have  not  achieved  a  conclusive  result  yet. 
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e.  Optimize  USE  code  for  using  the  human  data  (months  3-24).  I  have 
developed  a  novel  USE  technology  and  software  under  a  Breast  Cancer 
Research  Foundation  research.  The  work  thus  far  has  been  extremely 
promising:  a  manuscript  has  already  been  accepted  for  publication  in  the 
IEEE  Trans.  Med.  Imag.  I  will  be  further  improving  my  USE 
implementation  to  enhance  visualization  of  the  lumpectomy  bed  and 
ductal  tissue. 

Extremely  promising  results  have  been  obtained  in  this  area.  We  have 
developed  new  elastography  techniques  that  generate  superb  images  from  the 
human  data.  The  attached  MICCAI  paper  contains  some  of  the  details  of  the 
method.  Another  journal  paper  is  accepted  to  the  IEEE  TMI  (the  journal  with 
the  highest  impact  factor  in  medical  imaging).  And  finally,  a  paper  has  been 
submitted  to  IEEE  TMI  that  utilizes  multiple  images  to  calculate  strain 
images.  Since  the  paper  is  still  in  the  submission  stage,  we  cannot  attach  it  to 
this  report.  The  references  to  the  3  papers  are: 

Cl.  Rivaz,  H.,  Foroughi,  P.,  Fleming,  I.,  Zellars,  R.,  Boctor,  E.,  Hager,  G., 
Tracked  Regularized  Ultrasound  Elastography  for  Targeting  Breast 
Radiotherapy,  Medical  Image  Computing  and  Computer  Assisted 
Intervention,  MICCAI,  London,  UK,  Sept.  2009,  pp  507-515 

Jl.  Rivaz,  H.,  Boctor,  E.,  Choti,  M.,  Hager,  G.,  Real-Time  Regularized 
Ultrasound  Elastography,  IEEE  Trans.  Medical  Imaging  (Accepted,  in  press) 

J2.  Rivaz,  H.,  Boctor,  E.,  Choti,  M.,  Hager,  G.,  “Ultrasound  Elastography 
Using  Multiple  Images”,  IEEE  Trans.  Medical  Imaging  (submitted) 

Cl  and  Jl  are  attached.  J2  is  still  under  review. 

For  details  of  optimizing  the  USE  algorithm,  please  see  the  attachments. 


f.  Performance  optimization  and  refinements  of  subsystems  (months  18- 
36). 

The  ultrasound  elastography  subsystem  is  almost  finalized.  We  have  made  the 
code  also  available  to  the  public1,  so  that  groups  who  work  in  different 
applications  of  ultrasound  elastography  can  exploit  our  novel  and  high-quality 


1  Available  online  at  www.cs.ihu.edu/~rivaz/Ultrasound  Elastography 
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elastography  estimation  technique.  We  could  confidently  say  that  the  method 
can  be  implemented  commercially  with  small  modifications. 


2.  Obtain  tracked  US  and  USE  scans  of  the  patients  weekly  during  the 
radiotherapy  (months  2-14).  This  data  will  be  acquired  from  the  same  patients 
from  whom  US  and  USE  data  were  acquired  in  step  l.b. 

As  mentioned  in  the  lb  Section,  elastography  and  ultrasound  images  are  both 
acquired  from  patients. 


3.  Develop  US  to  US  registration  (month  25-27).  This  task  is  partly  solved  in  task 
l.c.  (US  to  CT  registration).  As  mentioned,  US  to  US  registration  is  one  of  the 
steps  required  in  US  to  CT  registration.  I  will  optimize  this  step  to  work  with 
breast  images,  as  opposed  to  simulation  images  in  task  l.c. 

As  mentioned  before,  registration  of  the  pseudo-ultrasound  to  the  ultrasound  images 
are  required  to  perform  CT  to  ultrasound  registration.  Please  refer  to  the  task  l.c  for 
more  details. 


KEY  RESEARCH  ACCOMPLISHMENTS: 

•  Development  of  a  novel,  real-time,  robust  and  high  quality  elastography 
techniques 

•  Development  of  a  tracked  ultrasound  system  using  magnetic  trackers 

•  Study  of  ultrasound  images  of  lumpectomy  cavity 

•  Study  of  elastography  images  of  lumpectomy  cavity 

•  Data  acquisition  from  more  than  15  patients 

•  Developed  the  2D  AM  method  which  calculates  high  quality  2D  strain  images  in 
real-time 

•  Introduced  ElastMI  (Elastography  using  Multiple  Images),  a  novel  method  which 
generates  high  quality  strain  images  by  utilizing  multiple  ultrasound  images. 


REPORTABLE  OUTCOMES: 

1.  Rivaz,  H.,  Foroughi,  P.,  Fleming,  I.,  Zellars,  R.,  Boctor,  E.,  Hager,  G.,  “Tracked 
Regularized  Ultrasound  Elastography  for  Targeting  Breast  Radiotherapy”,  Medical 
Image  Computing  and  Computer  Assisted  Intervention,  MICCAI,  London,  UK,  Sept. 
2009,  pp  507-515.  [Acceptance  rate:  34%]  [Awarded  MICCAI  Travel  Grant] 
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2.  Rivaz,  H.,  Boctor,  E.,  Choti,  M.,  Hager,  G.,  “Regularized  Ultrasound  Elastography”, 
IEEE  Trans.  Medical  Imaging  (highest  impact  factor  in  medical  imaging  engineering 
journals),  2011,  in  press. 

3.  Rivaz,  H.,  Boctor,  E.,  Choti,  M.,  Hager,  G.,  “Ultrasound  Elastography  Using  Multiple 
Images”,  IEEE  Trans.  Medical  Imaging  (submitted) 


Training: 

The  training  obtained  during  the  year  of  2010  was  mainly  in  the  following  forms: 

•  Reading  papers  related  to  partial  breast  radiation  therapy  techniques 

•  Helping  in  preparation  of  Institutional  Review  Boards  (IRB)  for  patient  trials 

•  Meetings  with  radiation  oncologists  on  the  problems  of  current  breast 
radiotherapy  workflow 

•  Analyzing  and  enhancing  the  CT,  ultrasound  and  strain  images  obtained  from  the 
lumpectomy 


CONCLUSION: 

Data  acquisition  from  lumpectomy  patients  went  smoothly.  We  acquired  the  ultrasound 
data  few  minutes  before  the  CT  scan  was  acquired.  This  means  that  we  did  not  alter  the 
current  medical  work-flow.  Ultrasound  data  acquisition  took  only  few  minutes, 
minimizing  the  cost  and  patient  discomfort.  Magnetic  tracking  provided  ease  of  use  and 
enough  reliability  for  registration  of  ultrasound  to  CT  images  and  for  reconstructing  3D 
volumes  from  2D  data.  The  novel  ultrasound  elastography  technique  has  many 
advantages  compared  to  the  previous  methods.  We  chose  the  novel  application  of  the 
lumpectomy  cavity  localization  as  the  hard  scar  tissue  around  the  lumpectomy  is 
relatively  thin  and  demands  a  high  resolution  elastography  method.  Also,  incoherent  fluid 
motions  in  the  cavity  causes  large  decorrelations,  requiring  a  robust  method.  The  merit  of 
adding  ultrasound  and  elastography  to  delineation  of  the  lumpectomy  bed  is  not  clear  to 
us  yet. 
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Real-Time  Regularized  Ultrasound  Elastography 

Hassan  Rivaz,  Emad  M.  Boctor,  Michael  A.  Choti  and  Gregory  D.  Hager 


Abstract — This  paper  introduces  two  real-time  elastography 
techniques  based  on  analytic  minimization  (AM)  of  regularized 
cost  functions.  The  first  method  (ID  AM)  produces  axial  strain 
and  integer  lateral  displacement,  while  the  second  method  (2D 
AM)  produces  both  axial  and  lateral  strains.  The  cost  functions 
incorporate  similarity  of  RF  data  intensity  and  displacement 
continuity,  making  both  AM  methods  robust  to  small  decorrela¬ 
tions  present  throughout  the  image.  We  also  exploit  techniques 
from  robust  statistics  to  make  the  methods  resistant  to  large 
local  decorrelations.  We  further  introduce  Kalman  filtering  for 
calculating  the  strain  field  from  the  displacement  field  given  by 
the  AM  methods.  Simulation  and  phantom  experiments  show 
that  both  methods  generate  strain  images  with  high  SNR,  CNR 
and  resolution.  Both  methods  work  for  strains  as  high  as  10% 
and  run  in  real-time.  We  also  present  in-vivo  patient  trials  of 
ablation  monitoring.  An  implementation  of  the  2D  AM  method 
as  well  as  phantom  and  clinical  RF-data  can  be  downloaded  from 
http://www.cs.jhu.edu/~rivaz/Ultrasound_Elastography/. 

Index  Terms — Real-time  ultrasound  elastography,  2D  strain, 
regulariation,  robust  estimation,  Kalman  filter,  RF  ablation. 


I.  Introduction 

ELASTOGRAPHY  involves  imaging  the  mechanical  prop¬ 
erties  of  tissue  and  has  numerous  clinical  applications. 
Among  many  variations  of  ultrasound  elastography  [1],  [2], 
[3],  [4],  our  work  focuses  on  real-time  static  elastography,  a 
well-known  technique  that  applies  quasi-static  compression  of 
tissue  and  simultaneously  images  it  with  ultrasound.  Within 
many  techniques  proposed  for  static  elastography,  we  focus 
on  freehand  palpation  elasticity  imaging  which  involves  de¬ 
forming  the  tissue  by  simply  pressing  the  ultrasound  probe 
against  it.  It  requires  no  extra  hardware,  provides  ease  of  use 
and  has  attracted  increasing  interest  in  recent  years  [5],  [6], 
[7],  [8],  [9],  [10].  Real-time  elastography  is  of  key  importance 
in  many  diagnosis  applications  [11],  [6],  [12],  [8],  [13]  and  in 
guidance/monitoring  of  surgical  operations  [14],  [15],  [16]. 

Global  and  local  decorrelation  between  the  pre-  and  post¬ 
compression  ultrasound  images  compromises  the  quality  of  the 
elasticity  images.  The  main  sources  of  global  decorrelation 
in  freehand  palpation  elastography  are  change  of  speckle 
appearance  due  to  scatterer  motion  and  out-of-plane  motion 
of  the  probe  (axial,  lateral  and  out-of-plane  directions  are 
specified  in  Figure  1).  Examples  of  local  decorrelation  are:  (1) 
a  decrease  in  the  ultrasonic  signal  to  noise  ratio  with  depth,  (2) 
low  correlation  close  to  arteries  due  to  complex  motion  and 
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inside  blood  vessels  due  to  blood  motion,  (3)  extremely  low 
correlation  in  lesions  that  contain  liquid  due  to  the  incoherent 
fluid  motion  [17],  [8],  and  (4)  out-of-plane  motion  of  movable 
structures  within  the  image  [17]. 

Most  elastography  techniques  estimate  local  displacements 
of  tissue  based  on  amplitude  correlation  [18],  [2]  or  phase 
correlation  of  the  radio-frequency  (RF)  echoes  [19],  [20], 
[21].  Assuming  a  stationary  signal  model  for  the  RF  data, 
the  use  of  large  correlation  windows  helps  to  reduce  jitter 
errors  (variance)  for  all  motion  field  estimation  techniques 
studied  in  [18],  [22].  This  is  intuitive  as  larger  windows 
contain  more  information.  However,  in  practice  RF  data  is 
not  stationary  and,  for  large  deformations,  the  decorrelation 
increases  with  window  size.  Therefore,  in  addition  to  reducing 
the  spatial  resolution  [23],  larger  windows  result  in  significant 
signal  decorrelation  [24],  [23],  [18].  Coarse-to-fine  hierarchi¬ 
cal  search  is  used  in  [23]  to  combine  the  accuracy  of  large 
windows  with  the  good  spatial  resolution  of  small  window. 
However,  the  issue  of  signal  decorrelation  within  the  window 
remains  unresolved  in  this  approach  and  can  cause  the  highest 
level  of  the  hierarchical  search  to  fail. 

All  of  the  aforementioned  methods  either  do  not  calculate 
the  lateral  displacement  or  they  just  calculate  an  approximate 
integer  lateral  displacement.  A  2D  displacement  field  is  re¬ 
quired  to  calculate  the  thermal  expansion,  lateral  and  shear 
strain  fields  [25]  (i.e.  reconstruct  the  strain  tensor),  Poisson’s 
ratio  and  Young’s  modulus  [26],  [27].  The  axial  resolutions  of 
ultrasound  is  determined  by  the  pulse  length,  and  the  lateral 
resolutions  is  dictated  by  the  center  frequency  of  the  excitation 
and  the  transducer  pitch.  Therefore,  the  lateral  resolution  is  of 
order  of  magnitude  lower  than  axial  resolution.  As  a  result, 
few  2D  elastography  techniques  have  been  proposed  to  date. 
Initially,  2D  motion  estimation  started  in  the  field  of  blood  flow 
estimation  using  speckle  tracking  [28].  Designed  for  blood 
flow  estimation,  these  techniques  are  not  immediately  suitable 
for  elastography  which  involves  tissue  deformation. 

Attaching  a  coordinate  system  to  the  ultrasound  probe  as 
in  Figure  1 ,  z,  x  and  y  in  the  ultrasound  image  are  generally 
defined  as  axial,  lateral  and  elevational  directions.  Assume  that 
the  applied  compression  to  the  tissue  is  the  Z  direction,  and 
attach  a  coordinate  system  X,  Y,  Z  to  the  applied  force.  Letting 
dz  and  djsr  be  the  displacements  in  the  Z  and  N  directions 
where  X_LZ,  axial  and  transverse  strains  are  ddz/dZ  and 
ddjy/dN.  In  most  experimental  setups  (including  freehand 
palpation  elastography),  z  and  Z  are  parallel  and  N  will 
be  either  lateral  or  out-of-plane,  and  therefore  d n  cannot  be 
estimated  accurately.  To  calculate  an  accurate  transverse  strain, 
Z  and  z  are  perpendicular  in  [29]  by  applying  the  compression 
force  perpendicular  to  the  ultrasound  imaging  axis.  Therefore, 
transverse  strain  is  in  the  z  direction  of  the  ultrasound  probe 
and  hence  can  be  measured  with  high  accuracy.  However, 
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such  an  experimental  setup  is  not  possible  in  many  medical 
applications.  Beam  steering  has  been  used  to  solve  this  issue 
[30].  In  freehand  palpation  elastography,  beam  steering  causes 
2  and  Z  to  be  unparallel,  so  that  a  component  of  the  dx  is 
in  the  2  direction.  The  steering  angle  determines  the  angle 
between  2  and  Z.  Unfortunately,  large  steering  angles  are 
required  to  obtain  accurate  estimates  of  lateral  strain,  which 
is  possible  in  phased  arrays  and  not  in  linear  arrays. 

Lateral  strain1  estimation  is  obtained  in  [31]  by  iteratively 
calculating  axial  strain,  companding  RF  data  and  interpolating 
in  the  lateral  direction.  In  another  work  [32],  tissue  defor¬ 
mation  is  modeled  by  locally  affine  transformations  to  obtain 
both  lateral  and  axial  strains.  Change  of  speckle  appearance  is 
taken  into  account  by  proposing  a  Lagrangian  speckle  model 
[33].  Although  they  provide  high  quality  lateral  strain,  these 
techniques  are  computationally  expensive  and  are  not  suitable 
for  real-time  implementation. 

Axial  strain  is  used  in  [34]  to  enhance  the  quality  of  lateral 
displacement  estimation.  Tissue  is  assumed  to  be  incompress¬ 
ible  and  isotropic  and  therefore  axial,  lateral  and  out-of-plane 
strains  should  add  to  zero.  However,  many  tissues  cannot  be 
considered  incompressible.  In  fact,  some  research  has  even 
focused  on  imaging  the  ratio  of  the  axial  and  lateral  strain 
(i.e.  the  Poisson’s  ratio  v)  [31]. 

While  most  previously  mentioned  methods  use  tissue  mo¬ 
tion  continuity  to  confine  the  search  range  for  the  neighboring 
windows,  the  displacement  of  each  window  is  calculated 
independently  and  hence  is  sensitive  to  signal  decorrelation. 
Since  data  alone  can  be  insufficient  due  to  signal  decorrelation, 
Pellot-Barakat  et  al.  [35]  proposed  minimizing  a  regularized 
energy  function  that  combines  constraints  of  conservation  of 
echo  amplitude  and  displacement  continuity.  In  another  work 

[36] ,  both  signal  shift  and  scale  are  found  through  minimiza¬ 
tion  of  a  regularized  cost  function.  The  computation  time  of 
these  methods  is  reported  to  be  few  minutes  and  therefore 
are  not  immediately  suitable  for  real  time  elastography.  In 

[37] ,  [38],  few  phase-based  methods  are  regularized  and  strain 
and  elasticity  modulus  images  are  obtained.  The  regulariza¬ 
tion  term  is  the  Laplacian  (second  derivative)  of  the  motion 
field  and  is  spatially  variant  based  on  the  peak-value  of  the 
correlation  coefficient.  The  regularization  makes  the  method 
significantly  more  robust  to  signal  decorrelation.  However,  it  is 
still  prone  to  decorrelation  within  each  window  especially  for 
large  strain  rates.  In  a  recent  work  [39],  a  displacement  field  is 
first  calculated  by  minimizing  phase  differences  in  correlation 
windows  [21].  The  strain  image  is  then  estimated  from  the 
displacement  field  by  optimizing  a  regularized  cost  function. 
The  regularization  assures  smooth  strain  image  calculation 
from  the  noisy  displacement  estimates. 

We  have  proposed  optimizing  a  recursive  regularized  cost 
function  using  Dynamic  Programming  (DP)  [40].  DP  is  used 
to  speed  the  optimization  procedure,  but  it  only  gives  in¬ 
teger  displacements.  Subsample  displacement  estimation  is 
possible  [40],  but  it  is  computationally  expensive,  particularly 
if  subsample  accuracy  is  needed  in  both  axial  and  lateral 

!We  hereafter  assume  the  applied  force  is  in  the  z  direction  (i.e.  Z  and  z 
are  parallel)  and  therefore  we  use  the  term  lateral  strain  instead  of  the  term 
transverse  strain. 


directions.  Therefore,  only  axial  subsample  displacement  is 
calculated  [40].  In  addition,  a  fixed  regularization  weight  is 
applied  throughout  the  image.  To  prevent  regions  with  high 
local  decorrelation  from  introducing  errors  in  the  displacement 
estimation  one  should  use  large  weights  for  the  regularization 
term,  resulting  in  over- smoothing. 

In  this  paper,  we  present  two  novel  real-time  elastogra¬ 
phy  methods  based  on  analytic  minimization  (AM)  of  cost 
functions  that  incorporate  similarity  of  echo  amplitudes  and 
displacement  continuity.  Similar  to  DP,  the  first  method  gives 
subsample  axial  and  integer  lateral  displacements.  The  second 
method  gives  subsample  2D  displacement  fields  and  2D  strain 
fields.  The  size  of  both  displacement  and  strain  fields  is  the 
same  size  as  the  RF-data  (i.e.  the  methods  are  not  window 
based  and  the  displacement  and  strain  fields  are  calculated 
for  all  individual  samples  of  RF-data).  We  introduce  a  novel 
regularization  term  and  demonstrate  that  it  minimizes  displace¬ 
ment  underestimation  caused  by  smoothness  constraint.  We 
also  introduce  the  use  of  robust  statistics  implemented  via 
iterated  reweighted  least  squares  (IRLS)  to  treat  uncorrelated 
ultrasound  data  as  outliers.  Finally,  for  the  first  time  to  the  best 
of  our  knowledge  we  introduce  the  use  of  Kalman  filtering 
[41]  for  calculating  strain  image  from  the  displacement  field. 
Simulation  and  experimental  results  are  provided  for  quan¬ 
titative  validation.  The  paper  concludes  with  a  clinical  pilot 
study  utilizing  this  system  for  monitoring  thermal  ablation  in 
patients  with  liver  tumors. 

II.  Methods 

Assume  that  the  tissue  undergoes  a  deformation  and  let  I\ 
and  I 2  be  two  images  acquired  from  the  tissue  before  and  after 
the  deformation.  Letting  I\  and  I2  be  of  size  m  x  n  (Figure 
1),  our  goal  is  to  find  two  matrices  A  and  L  where  the  (i,j)th 
component  of  A  (a^j)  and  L  ( lij )  are  the  axial  and  lateral 
motion  of  the  pixel  (i,j)  of  Ii  (we  are  not  calculating  the 
out-of-plane  motion).  The  axial  and  lateral  strains  are  easily 
calculated  by  spatially  differentiating  A  in  the  axial  direction 
(resulting  in  Aa)  and  L  in  the  lateral  direction  (resulting  in 
L[).  The  shear  strains  (not  calculated  in  this  work)  can  also 
be  easily  calculated  by  spatially  differentiating  A  in  the  lateral 
direction  (resulting  in  Ai)  or  L  in  the  axial  direction  (resulting 
in  La). 

In  this  section,  we  first  give  a  brief  overview  of  a  previous 
work  (DP)  that  calculates  integer  values  for  A  and  L.  We  then 
propose  ID  Analytic  Minimization  (AM)  as  a  method  that 
takes  the  integer  displacement  field  from  DP  and  refines  the 
axial  displacement  component.  We  then  introduce  2D  Analytic 
Minimization  (AM)  that  takes  the  integer  displacement  of  a 
single  RF-line  from  DP  and  gives  the  subsample  axial  and 
lateral  displacement  fields  for  the  entire  image.  We  conclude 
this  section  by  presenting  a  technique  for  calculating  smooth 
strain  field  from  the  displacement  field  using  Kalman  filtering. 

A.  Dynamic  Programming  (DP) 

In  order  to  present  the  general  DP  formulation,  we  consider 
a  single  column  j  (an  RF-line)  in  I\  (the  image  before 
deformation)  in  Figure  1.  Let  m  and  n  be  the  length  of  the 
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z  (axial) 


- ►  /7,  before  deformation 

1  2  -  -  -  j  j+1  -  -  -  n 


G,j) 

- ►  I2,  after  deformation 

1  2  *  *  *  j  j+1  •  •  •  n 


( 

X 

[+\j  ’  j+1i,j) 

1 


Fig.  1.  Axial,  lateral  and  out-of-plane  directions.  The  coordinate  system  is  attached  to  the  ultrasound  probe.  The  sample  (i,j)  marked  by  x  moved  by 
(aij ,  hj).  aij  and  lij  are  respectively  axial  and  lateral  displacements  and  initially  are  integer  in  DP. 


RF-lines  and  the  number  of  RF-lines  in  the  images  (Figure  1). 
Let  di  and  U  denote  the  axial  and  lateral  displacements  of  the 
ith  sample  of  the  RF-line  in  column  j.  In  DP  elastography 
[40],  a  regularized  cost  function  is  generated  by  adding  the 
prior  of  displacement  continuity  (the  regularization  term)  to 
an  amplitude  similarity  term.  The  displacement  continuity  term 
for  column  j  is 

Rj  i&i  i  li  i  — 1 •>  l)  —  *Ta  ((2^  Q>i— l)  +  (Ji  li— 1 )  (1) 


which  forces  the  displacements  of  the  sample  i  (i.e.  ai  and  /*) 
be  similar  to  the  displacements  of  the  previous  sample  i— 1  (i.e. 
di-i  and  k- 1).  a a  and  ai  are  axial  and  lateral  regularization 
weights  respectively.  We  write  Rj{ai ,  li ,  a*_i,  k-±)  to  indicate 
the  dependency  of  di  and  li  on  j.  The  regularized  cost  function 
for  column  j  is  then  generated  as  following 


Cj(di,  )  —  [Ii(i,j)  —  I2 (J*  +  cii,j  +  k)]2  + 


da,di 


Cj(da,  di,i  1)  +  Cj_ 1  {da ,  di ,  i) 


RA 


j  i  ^5  da  5  di )  ^ 


(2) 


where  da  and  di  are  temporary  displacements  in  the  axial  and 
lateral  directions  that  are  varied  to  minimize  the  term  in  the 
bracket.  After  calculating  Cj  for  i  =  2  •  •  •  m,  Cj  is  minimized 
at  i  —  m  giving  am  and  Zm.  The  and  li  values  that  have 
minimized  the  cost  function  at  i  =  m  are  then  traced  back 
to  i  =«  1,  giving  integer  ai  and  li  for  all  samples  of  jth  line. 
The  process  is  performed  for  the  next  line  j  +  1  until  the 
displacement  of  the  whole  image  is  calculated.  The  2D  DP 
method  gives  integer  axial  and  lateral  displacement  maps.  In 
[40],  we  performed  hierarchical  search  to  obtain  subsample 
axial  displacement  (the  lateral  displacement  was  not  refined  to 
subsample).  DP  is  an  efficient  method  for  global  optimization 
and  has  been  used  extensively  in  many  applications  in  com¬ 
puter  vision  including  solving  for  optimal  deformable  models 
[42] .  In  the  next  section,  we  propose  an  alternative  method  for 
calculating  subsample  axial  displacement  which  is  both  faster 
and  more  robust  than  hierarchical  DP. 


B.  ID  Analytic  Minimization  (AM) 

Tissue  deformations  in  ultrasound  elastography  are  usually 
very  small  and  therefore  a  subsample  displacement  estima¬ 


tion  is  required.  We  now  develop  a  method  that  analytically 
minimizes  a  regularized  cost  function  and  gives  the  refined 
displacement  field  following  the  work  presented  in  [16].  We 
first  consider  a  specialization  of  Equation  2  in  which  we  only 
consider  refining  axial  displacements  to  subsample  level. 

Having  the  integer  displacements  a*  and  li  from  DP,  it 
is  desired  to  find  A  ai  values  such  that  ai  +  A  gives  the 
value  of  the  displacement  at  the  sample  i  for  i  =  1  •  •  •  m 
(li,  ai  and  A  a*  correspond  to  line  j.  Hereafter,  wherever 
the  displacements  correspond  to  the  jth  line,  j  is  omitted  to 
prevent  notation  clutter).  Such  A values  will  minimize  the 
following  regularized  cost  function 


Ci(Aa1,-..,Aam)  =  E£1{ 

[h (hj)  ~  +  Q>i  +  Aai,  j  +  li)f  + 

+  A&2  Q>i—  1  Aa$_i)  + 

OLi(di  +  A  ai  -  dij- 1  -  Aaij-1)2]}  (3) 

where  aa  >  0  and  ai  >  0  are  tunable  axial  and  lateral 
regularization  weights  and  subscript  j  —  1  refers  to  the  previous 
RF-line  (adjacent  RF-line  in  the  lateral  direction).  Substituting 
l^(i  +  di  +  Adi)  with  its  first  order  Taylor  expansion  approx¬ 
imation  around  di,  we  have 

Ci(Aa1,...,Aam)  =  E£1{ 

\Ii(hj)  ~  I2 (4  +  di,  j  +  li)  —  I'2(i  +  cii,  j  +  li)Adi)]  + 

aa(ai  +  Aa^  —  ai- 1  —  Aai-i)2  + 

oti(ai  +  A  ai  -  dij- 1  -  Aa^y-i)2]}  (4) 

where  I2  is  the  derivative  of  the  I2  in  the  axial  direction.  The 
optimal  Aai  values  occur  when  the  partial  derivative  of  Cj 
with  respect  to  Aai  is  zero.  Setting  =  0  for  i  =t  1  •  •  •  m 
we  have 

(I2  +^0+^1) Aay  =  V2e-(aa'D-\-aii)sLj-\-aiSLj-U  (5) 
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[0  0  •••  0  -1  1  J 

where  I2  =  diag(V2(\  +  d\,j  +  k)  •  •  •  I^m  +  dm,j  +  k))9 
=  [Aoq  j  •  •  •  Aum?j]  ,  e  =  [ei  •  •  •  em]  ,  =  /i(i,  j) 

12(1  + di,j  +  h),  a  j  =  [aij  •  •  •  amj]T,  i  is  the  identity  matrix 
and  aj_i  is  the  total  displacement  of  the  previous  line  (i.e. 
when  the  displacement  of  the  j  —  1th  line  was  being  calculated, 
8Lj~i  was  updated  with  3-j-i  +  Aaj_i).  I2,  D  and  I  are 
matrices  of  size  m  x  m  and  Aa,  e  and  a  are  vectors  of  size 


Comparing  ID  AM  (as  formulated  in  Equation  5)  and  2D 
DP,  they  both  optimize  the  same  cost  function.  Therefore,  they 
give  the  same  displacement  fields  (up  to  the  refinement  level  of 
the  DP).  In  the  next  two  subsections,  we  will  further  improve 
ID  AM. 

1)  Biasing  the  Regularization :  The  regularization  term 

Oia (cii  +  Adi  —  di- 1  —  Aai-i)2  penalizes  the  difference 
between  a*  +  A  a*  and  a^_  1  +  Aa^_i,  and  therefore  can 
result  in  underestimation  of  the  displacement  field.  Such 
underestimation  can  be  prevented  by  biasing  the  regulariza¬ 
tion  by  e  to  aa(ai  +  A  a*  —  a^_  1  —  Aa^_i  —  e)2,  where 
e  =  (am  —  ai)/(m  —  1)  is  the  average  displacement  difference 
(i.e.  average  strain)  between  samples  i  and  i  —  1.  An  accurate 
enough  estimate  of  —  d\  is  known  from  the  previous 
line.  With  the  bias  term,  the  R.H.S.  of  Equation  5  becomes 
I2e  —  (cyaD  +  aiTjdLj  +  +  Aa^-i)  +  b  where  the  bias 

term  is  b  =  OLa[— e  0---0  e]T  (only  the  first  and  the  last 
terms  are  nonzero)  and  all  other  terms  are  as  before.  In  the 
other  words,  except  for  the  first  and  the  last  equations  in  this 
system,  all  other  m  —  2  equations  are  same  as  Equation  5. 

Equation  5  can  be  solved  for  Aa  j  in  4m  operations  since 
the  coefficient  matrix  r22  +  a  aD  +  rfyi  is  tridiagonal.  Utilizing 
its  symmetry,  the  number  of  operations  can  be  reduced  to  2m. 
The  number  of  operations  required  for  solving  a  system  with 
a  full  coefficient  matrix  is  more  than  m3/ 3,  significantly  more 
than  2m. 

2)  Making  Elastography  Resistant  to  Outliers:  Even  with 
pure  axial  compression,  some  regions  of  the  image  may  move 
out  of  the  imaging  plane  and  decrease  the  decorrelation.  In 
such  parts  the  weight  of  the  data  term  in  the  cost  function 
should  be  reduced.  The  data  from  these  parts  can  be  regarded 
as  outliers  and  therefore  a  robust  estimation  technique  can 
limit  their  effect.  Before  deriving  a  robust  estimator  for  Aa, 
we  rewrite  Equation  4  as 


C(Aa)  =  S-1p(ri)  +  i?(Aa) 


(7) 


where  =  I\{i)  —  hii  +  di)  —  I'2(i  +  a^)Aa^  is  the  residual, 
p(ri)  =  r\  and  R  is  the  regularization  term.The  M-estimate 
of  Aa  is  Aa  =  arg  minAa  {H  •A1p(n)  +  R( Aa)}  where  p(ri) 
is  a  robust  loss  function  [43].  The  minimization  is  solved  by 


setting  &  = 0: 


An) 


dr 


+  dR(  Aa)  =  Q 


dA  a*  dAdj 


(8) 


A  common  next  step  [44]  is  to  introduce  a  weight  function 
w,  where  w(ri).ri  =  This  leads  to  a  process  known 

as  “iteratively  reweighted  least  squares”  (IRLS)  [45],  which 
alternates  steps  of  calculating  weights  w{ri)  for  =  1  •  •  •  m 
using  the  current  estimate  of  Aa  and  solving  Equation  8  to 
estimate  a  new  Aa  with  the  weights  fixed.  Among  many 
proposed  shapes  for  w(-),  we  compared  the  performance  of 
Huber  [44],  [43] 


w 


1 

T 


N  <  t 

\n\>T 


and  Cauchy  [45] 


w(ri) 


1 

1  +  (U/T)2 


(9) 

(10) 


functions  and  discovered  that  the  more  strict  Cauchy  function 
(which  decreases  with  inverse  of  the  square  of  the  residual) 
is  more  suitable  in  our  application.  To  better  discriminate 
outliers,  we  calculate  the  residuals  at  linear  interpolation 
of  the  integer  sample  displacements  provided  by  DR  With  the 
addition  of  the  weight  function,  Equation  8  becomes 

(wI22+aaD+rfyi)Aaj  =  wl2e— (o;aD-ho;/i)a^+a/aj_^+b 

(11) 

where  w  =  diag{w(r\)  •  •  •  w(rm)).  This  equation  will  con¬ 
verge  to  a  unique  local  minimum  after  few  iterations  [45]. 
The  convergence  speed  however  depends  on  the  choice  of 
T,  which  in  this  work  is  defined  manually.  Since  the  Taylor 
approximation  gives  a  local  quadratic  approximation  of  the 
original  non-quadratic  cost  function,  the  effect  of  higher  orders 
terms  increase  if  A dj  is  large.  Assuming  that  DP  gives  the 
correct  displacements,  ||Aaj||  <  e  where  H-H^  is  the  infinity 

norm  and  e  <  0.5.  In  practice,  however,  e  <<  0.5  because  the 
linear  interpolation  of  the  DP  displacements  (which  is  very 
close  to  the  correct  displacement)  is  used  to  calculate  the 
residuals  r^.  Therefore,  a  small  value  can  be  assigned  to  T 
in  ID  AM  provided  that  DP  results  are  trusted. 

The  coefficient  matrix  Q  =  wl22  +  aaE>  +  ail  in  Equation 
11  is  the  Hessian  of  the  cost  function  C  whose  minimum 
is  sought.  This  matrix  is  strictly  diagonally  dominant  (i.e. 
\qu\  >  \qij |  for  all  i  where  qij  is  the  z,  jth  element  of 
Q ),  symmetric  and  all  diagonal  entries  are  positive.  Therefore, 
it  is  positive  definite,  which  means  that  setting  the  gradient 
of  C  to  zero  results  in  the  global  minimum  of  C  (not  in  a 
saddle  point,  a  local  maximum  or  a  local  minimum).  All  of 
the  ID  AM  results  presented  in  this  work  are  obtained  with 
one  iteration  of  the  above  equation. 

ID  AM  takes  the  integer  axial  and  lateral  displacement 
fields  from  DP  and  gives  refined  axial  displacement.  It  inherits 
the  robustness  of  DP  and  adds  more  robustness  when  calculat¬ 
ing  the  fine  axial  displacements  via  IRLS.  However,  there  are 
redundant  calculations  in  this  method  which  are  eliminated  in 
2D  AM  as  described  next. 


C.  2D  Analytic  Minimization  (AM) 

In  2D  AM,  we  modify  Equation  2  to  calculate  subsample 
axial  and  lateral  displacement  fields  simultaneously.  The 
outline  of  our  proposed  algorithm  is  as  follows 
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1)  Calculate  the  integer  axial  and  lateral  displacements  of 
one  or  more  seed  RF-lines  (preferably  in  the  middle  of 
the  image)  using  DP  (Equation  2).  Calculate  the  linear 
interpolation  of  the  integer  displacements  as  an  initial 
subsample  estimate. 

2)  Calculate  subsample  axial  and  lateral  displacements  of 
the  seed  RF-line  using  2D  AM,  as  explained  below. 
Add  the  subsample  axial  and  lateral  displacements  to 
the  initial  estimate  to  get  the  displacement  of  the  seed 
line. 

3)  Propagate  the  solution  to  the  right  and  left  of  the  seed 
RF-line  using  the  2D  AM  method,  taking  the  displace¬ 
ment  of  the  previous  line  as  the  initial  displacement 
estimate. 

Benefits  of  2D  AM  are  two-fold.  First  it  computes  subsam¬ 
ple  displacements  in  both  axial  and  lateral  directions.  Lateral 
strain  contains  important  information  from  tissue  structure  that 
is  not  available  from  axial  strain  [31],  [46],  [47].  Second,  it 
is  only  required  to  calculate  the  displacement  of  a  single  line 
using  DP  (the  seed),  eliminating  the  need  to  have  the  integer 
displacement  map  for  the  entire  image.  This  is  significant  as  in 
the  ID  AM  method,  the  initial  step  to  calculate  the  2D  integer 
displacements  using  DP  takes  about  10  times  more  than  the 
ID  AM. 

Assume  that  initial  displacement  estimates  in  the  axial 
direction,  a$,  and  in  the  lateral  direction,  k,  are  known  for 
all  i  =  1  •  •  •  m  samples  of  an  RF-line.  Note  that  a*  and  U  are 
not  integer;  for  the  seed  line  they  are  the  linear  interpolation 
of  the  integer  DP  displacements  and  for  the  rest  of  the  lines 
are  the  displacement  of  the  previous  line.  It  is  desired  to  find 
Acti  and  A li  values  such  that  the  duple  (a*  +  A ai,  U  +  A li) 
gives  the  axial  and  lateral  displacements  at  the  sample  i.  Such 
(A di,  Aa^  values  will  minimize  the  following  regularized  cost 
function 


Cj ( Atti ,  •  •  • ,  Aam  ,  All,  *  *  *  ^  = 

T^=i{[h{i,j)  —  I2(i  +  di  +  A  di,j  +  U  +  A  li)]  + 

a(di  +  Aa^  —  a^_  i  —  Aa^_i)2  + 

(da  (Ji  +  A  li  —  li- 1  —  Ali-i)2  +  Pi(li  +  A  li  —  hj— i)2}  (12) 

where  I(i,  j)  is  the  ith  sample  on  the  jth  RF-line.  Since  we 
perform  the  calculations  for  one  RF-line  at  a  time,  we  dropped 
the  index  j  to  simplify  the  notations:  di,  li ,  Adi  and  A  li  are 
dij,  lij ,  A dij  and  A lij.  kj-i  is  the  lateral  displacement 
of  the  previous  RF-line  (note  that  hj-i  is  the  total  lateral 
displacement  of  the  previous  line,  i.e.  when  the  displacement 
of  the  j  —  1th  line  was  being  calculated,  kj-i  was  updated 
with  lij-i~\-  Alij-i).  Since  in  the  first  iteration  a*  and  li  (the 
initial  displacement  estimates)  are  in  fact  the  displacements  of 
the  previous  RF-line,  for  the  first  iteration  we  have  kj~i  =  h- 
This  simplifies  the  last  term  in  the  R.H.S.  to  p[A l2.  The 
regularization  terms  are  a ,  (3a  and  (3\\  a  determines  how 
close  the  axial  displacement  of  each  sample  should  be  to 
its  neighbor  on  the  top  and  pa  and  (3[  determine  how  close 
lateral  displacement  of  each  sample  should  be  to  its  neighbors 
on  the  top  and  left  (or  right  if  propagating  to  the  left). 
If  the  displacement  of  the  previous  line  is  not  accurate,  it 
will  affect  the  displacement  of  the  next  line  through  the  last 


term  in  the  R.H.S.  of  Equation  12.  Although  its  effect  will 
decrease  exponentially  with  j,  it  will  propagate  for  few  RF- 
lines.  Therefore  we  set 


Pi 


Pi 

1  +  \rij- 1| 


(13) 


to  prevent  such  propagation  where  i  is  the  residual 
associated  with  the  displacement  of  the  ith  sample  of  the 
previous  line.  A  large  residual  indicates  that  the  displacement 
is  not  accurate  and  therefore  its  influence  on  the  next  line 
should  be  small,  which  is  realized  via  the  small  weight  (3[.  This 
is,  in  principle,  similar  to  guiding  the  displacement  estimation 
based  on  a  data  quality  indicator  [48].  The  effect  of  the  tunable 
parameters  a ,  pa  and  Pi  is  studied  in  the  Results  section. 
Writing  the  2D  Taylor  expansion  of  the  data  term  in  Equation 
12  around  ( i  +  di,j  +  li): 


1 2  ('l  +  di  +  A  di,j  +  li  %  A  li)  ~ 

I2  (p  +  Ui,j  +  li)  +  Adilya  +  Alil'2  l  (14) 

where  I2  a  and  I2  t  are  the  derivatives  of  the  I2  at  point 
(i  +  di,  j  +  li)  in  the  axial  and  lateral  directions  respectively. 
Note  that  since  the  point  (i  +  di,j  +  li)  is  not  on  the 
grid  {di  and  U  are  not  integer),  interpolation  is  required  to 
calculate  I2  a  and  I2  v  We  propose  a  method  in  Section  II-C1 
that  eliminates  the  need  for  interpolation.  The  optimal  (A 

A  li)  values  occur  when  the  partial  derivatives  of  Cj  with 

dC  ■  J 

respect  to  both  Adi  and  A U  are  zero.  Setting  -7^-  =  0  and 

dC  ■  % 

=  0  for  i  =  1  •  •  •  m  and  stacking  the  2m  unknowns  in 


Ad  =  [Aai  Ah 

Aa2  AZ2  •  • 

Aam 

A  lm 

T  and  the  2m  initial 

estimates  in  d  = 

[ai  h 

d2  l2 

•  drn 

Im]7 

we  have 

d'2 

CVi 

+  V2)Ad  = 

X2'e 

-vxd, 

(15) 
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where  =  didg( 0,  P[,  0,  P[,  •  •  • ,  0,  p[)  is  a  diagonal  matrix 
of  size  2m  x  2m  ,  X2  =  didg(3'2(l)  •  •  •  3'2(m))  is  a 
symmetric  tridiagonal  matrix  of  size  2m  x  2m  with 


y\i) 


V  2 

J2  ,a 

T'  T' 

12,a12,l 


T'  T' 

12,a12,l 

V  2 

i2 ,1 


(16) 


blocks  on  its  diagonal  entries  where  I2  a  and  I2  t  are  the 
derivatives  of  the  I2  at  point  (i  +  di,  j  +  li)  in  the  axial  and 
lateral  directions, 


X ^  =  dia5(40(l),  7^,(1),  40(2),/^(2).../'i0(m),/^(m)) 

(17) 

where  I2  a{i)  and  I2,i(i)f  are  calculated  at  point  (i  +  a^,  j  +  k), 
and  e  =  [ex  ex  e2  e2  •  •  -em]T,  e»  =  h(i,  j)  - I2(i  +  di,  j +  k). 

We  make  four  modifications  to  Equation  15:  First,  we  take 
into  account  the  attenuation  of  the  ultrasound  signal  with 
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depth.  As  the  signal  gets  weaker  with  depth,  the  first  term 
in  the  R.H.S.  of  Equation  15  (X2c)  gets  smaller.  This  results 
in  increasing  the  share  of  the  regularization  term  in  the  cost 
Cj  and  therefore  over- smoothing  the  bottom  of  the  image. 
The  attenuation  of  the  ultrasound  signal  [49]  reflected  from  the 
depth  d  is  ((d)  =  e-21°s(10)at/o^/20  wjiere  at  [ s  frequency 
dependent  attenuation  coefficient  of  tissue  and  is  equal  to  0.63 
dB/cm/MHz  for  fat  [49],  /o  is  the  center  frequency  of  the  wave 
(in  MHz)  and  d  is  in  cm.  Having  the  exponential  attenuation 
equation,  the  attenuation  level  at  sample  i  will  be 

1540xl02at/0  log(10) 

=  x  =  e  20/  s  x  io®  5  i  =  l---m  (18) 

where  1540  x  102  is  the  speed  of  sound  in  tissue  (in  cm/sec) 
and  fs  is  the  sampling  rate  of  the  ultrasound  system  (in 
MHz).  This  is  assuming  that  the  TGC  (time  gain  control)  is 
turned  off.  Otherwise,  the  TGC  values  should  be  taken  into 
account  in  this  equation.  Let  the  2 m  x  2 m  diagonal  matrix 
Z  be  Z  =  diag(Ci,Ci,C2,C2---Cm,Cm)-  To  compensate  for 
the  attenuation,  we  multiply  the  V\  and  V 2  matrices  in 
Equation  15  by  Z,  and  therefore  reduce  the  regularization 
weight  with  depth.  As  we  will  show  in  Sections  III  and 
IV,  the  regularization  weight  can  vary  substantially  with  no 
performance  degradation.  Therefore  approximate  values  of  the 
speed  of  sound  and  attenuation  coefficient  will  suffice.  Second, 
we  add  a  bias  term  in  the  regularization  similar  to  the  ID 
case.  Here  we  only  bias  the  axial  displacement  since  the 
difference  between  the  lateral  displacements  of  the  points  on 
a  RF-line  is  very  small,  usually  less  than  4  RF-lines.  Third, 
we  exploit  the  fact  that,  because  the  tissue  is  in  contact  with 
the  ultrasound  probe,  the  axial  displacement  of  the  top  of  the 
image  is  zero  relative  to  the  probe  (the  lateral  displacement  of 
the  top  of  the  image  is  not  zero  as  tissue  might  slip  under  the 
probe).  Therefore,  we  enforce  the  axial  displacement  of  the 
first  sample  to  be  zero  by  changing  the  first  row  of  X>i,  X2 
and  I2  .  Fourth,  we  make  the  displacement  estimation  robust 
via  IRLS  using  the  Cauchy  function  (Equation  10).  Similar 
to  ID  AM,  T  is  selected  manually.  For  the  first  (seed)  RF- 
line,  a  small  value  can  be  selected  for  T  if  DP  results  are 
trusted.  For  the  next  lines,  the  value  of  Ad  determines  the 
accuracy  of  the  Taylor  expansion  14:  for  a  small  Ad,  the 
residuals  of  the  inliers  are  small  and  therefore  a  small  T  can 
be  chosen,  while  for  a  large  Ad  the  inliers  might  give  large 
residuals  and  therefore  a  large  value  for  T  is  required.  Since 
the  tissue  motion  is  mostly  continuous,  Ad  mostly  depends 
on  the  lateral  sampling  of  the  image  (i.e.  the  number  of  A-line 
per  cm).  Therefore  if  many  A-lines  are  given  per  cm  of  the 
image  width,  a  small  value  of  T  will  give  the  optimum  results. 
Since  the  amplitude  of  signal  is  decreasing  due  to  attenuation, 
we  decrease  the  IRLS  parameter  T  with  depth  by  multiplying 
it  with  Q  at  each  sample  i.  With  these  modifications,  Equation 
15  becomes 

(yVl'2  +  ZDi  +  ZD2)Ad  =  Wife  -  ZDjd  +  s  (19) 

where  W  =  diag(0,  w(n),  w(r2),  w(r2)  •  •  •  w(rm),  w(rm)) 
(i.e.  W2j)2j  =  W2j— ij2j— i  =  w(ri)  for  i  =  1  •  •  •  m  except  for 
Wi,i  =  0  which  guarantees  the  displacement  of  the  first  sam¬ 
ple  to  be  zero)  is  the  weight  function  determined  by  the  resid- 
uals  r*  =  h(i,j)  -  [l2(i  +  dt,j  +  af  +  A djfz  +  Aajfx], 


w  is  as  before  (Equation  10),  the  bias  term  s  is  a  vector  of 
length  2m  whose  all  elements  are  zero  except  the  2m  —  1th 
element: 82m- 1  =  ae,  and  e  =  (drn  —  d\)/(m  —  1)  is 
as  before.  Similar  to  Equation  11,  the  coefficient  matrix 
Q  =  WX2  +  ZD/  +  ZV2  is  strictly  diagonally  dominant, 
symmetric  and  all  the  diagonal  entries  are  positive.  Therefore 
Q  is  positive  definite  which  means  that  solving  Equation  19 
results  in  the  global  minimum  of  the  cost  function  C.  The 
updated  displacement  field  (axial  and  lateral)  will  be  d  +  Ad. 

Equation  19  can  be  solved  for  Ad  in  9m  operations  since 
the  coefficient  matrix  WX2  +  ZT*i  +  ZD2  is  pentadiagonal 
and  symmetric.  This  number  is  again  significantly  less  than 
(2m)3/3,  the  number  of  operations  required  to  solve  a  full 
system. 

1 )  Inverse  Gradient  Estimation:  With  the  subsample  initial 
displacement  field,  the  Taylor  expansion  should  be  written 
around  off-grid  points,  which  requires  calculation  of  image 
gradient  at  these  points  (matrices  X2  and  X2  in  Equation  19). 
In  Figure  2  (a),  this  is  equivalent  to  calculating  gradient  of  /2 
on  the  off-grid  marks.  There  are  two  disadvantages  associated 
with  this:  1)  it  requires  interpolation  of  the  gradients,  and  2) 
the  image  gradient  should  be  recalculated  after  each  iteration. 
As  proposed  by  [44],  [50],  image  gradient  can  be  instead 
calculated  at  on-grid  locations  on  image  1  in  the  following 
way. 

Consider  two  problems:  (1)  to  find  the  matches  for  grid 
points  on  //  having  the  initial  off-grid  estimates  on  /2,  and 
(2)  to  find  the  matches  for  the  off-grid  points  on  /2  having  the 
initial  grid  estimates  on  I\.  For  both  problems,  /2  values  must 
be  interpolated  on  the  off-grid  values.  However,  the  second 
problem  does  not  require  interpolation  of  the  image  gradient 
since  the  Taylor  expansion  is  written  around  grid  points  of 
1 1  (Figure  2  (b)).  It  is  shown  in  [51]  that  the  two  techniques 
converge  to  the  same  results.  Therefore,  on  one  hand  inverse 
gradient  calculation  is  both  faster  and  easier  to  implement, 
and  on  the  other  hand  it  causes  no  performance  degradation. 
Exploiting  this,  Equation  19  becomes 

(Wlf  +  Z£>1  +  ZD2)Ad  =  yyi/e  -  ZDid  +  s  (20) 

where  if  and  1[  are  now  calculated  on  the  grid  points  of 
image  1. 

All  the  2D  AM  results  presented  in  this  work  are  obtained 
using  Equation  20.  For  the  seed  line  where  the  initial  estimate 
might  be  inaccurate,  this  equation  is  iterated  multiple  times 
(about  10  times).  For  all  other  lines  this  equation  is  iterated 
only  once. 

D.  Strain  Estimation  Using  Kalman  Filter 

Strain  estimation  requires  spatial  derivation  of  the  displace¬ 
ment  field.  Since  differentiation  amplifies  the  signal  noise, 
least  squares  regression  techniques  are  commonly  used  to 
obtain  the  strain  field.  Adjacent  RF-lines  are  usually  processed 
independently  in  strain  calculation.  However,  the  strain  value 
of  each  pixel  is  not  independent  from  the  strain  value  of  its 
neighboring  pixels.  The  only  exception  is  the  boundary  of  two 
tissue  types  with  different  mechanical  properties  where  the 
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(a)  Three  samples  on  I\  (left)  and  corresponding  matches  on  1 2  (right)  (b)  Inverse  gradient  estimation 

Fig.  2.  Schematic  plot  of  subsample  tissue  motion,  (a)  In  I2  the  initial  estimates  (in  black)  are  updated  by  the  arrows  (three  components  of  Ad)  to  new 
estimates  (in  red)  after  an  iteration  of  2D  AM.  To  find  Ad  using  Equation  19,  it  is  required  to  calculate  image  gradient  at  the  off-grid  initial  estimate  locations 
(in  black)  on  I2.  (b)  Schematic  plot  of  two  RF-data  I\  and  I2,  each  sampled  at  three  locations  (black  dots).  The  black  dashed-dotted  arrow  shows  A  a  of 
the  sample  on  I\  (ignoring  the  regularization  term)  which  requires  calculating  the  gradient  on  I2  at  an  off-grid  location.  The  blue  dashed  arrow  shows  A  a 
of  an  off-grid  sample  on  I2  (ignoring  the  regularization  term)  which  requires  calculating  the  gradient  on  I\  at  an  on-grid  location.  Ignoring  second  order 
derivatives,  the  length  of  the  two  arrows  is  equal. 


strain  field  is  discontinuous.  We  use  the  prior  of  piecewise 
strain  continuity  via  a  Kalman  filter  to  improve  the  quality 
of  strain  estimation.  Although  locations  with  strain  disconti¬ 
nuity  are  limited,  we  will  develop  a  technique  to  take  such 
discontinuities  into  account. 

We  first  calculate  the  strain  using  least  squares  regression. 
Each  RF-line  is  first  differentiated  independently:  for  each 
sample  i,  a  line  is  fitted  to  the  displacement  estimates  in  a 
window  of  length  2k  +  1  around  i,  i.e.  to  the  samples  i  —  k 
to  i  +  k.  The  slope  of  the  line,  ,  is  calculated  as  the 
strain  measurement  at  i.  The  center  of  the  window  is  then 
moved  to  i  +  1  and  the  strain  value  Zi+ij  is  calculated.  We 
reuse  overlapping  terms  in  calculation  of  Zij  and  Zi+ ij,  and 
therefore  the  running  time  is  independent  of  the  window  length 
2fc+l.  Having  Zij  for  i  =  1  •  •  •  m  and  j  —  1  •  •  •  n,  we  propose 
the  following  algorithm  based  on  Kalman  filter  to  take  into 
account  the  prior  of  strain  continuity. 

Zij  are  the  noisy  measurements  of  the  underlying  strain  field 
6ij.  Since  the  z^j  values  are  calculated  using  axial  windows, 
we  apply  the  Kalman  filter  in  the  lateral  direction.  Let  rj  be  the 
Gaussian  process  noise  and  Sj  be  the  Gaussian  measurement 
noise  to  be  removed.  We  have  [52],  [41] 


ei,j  —  €i,J- 1  +  ri,j  (21) 

—  G,j  T~  &i,j  (22) 

Let  (note  the  super  minus)  be  our  a  priori  strain  estimate 
from  the  process  prior  to  step  j  (i.e.  from  the  Equation  21) 
and  iij  be  our  a  posteriori  strain  estimate  at  step  j  given 
measurement  Zj.  Let  also  the  variances  of  e^.  and  be 
respectively  p~  and  p.  The  time  update  (i.e.  prior  estimation) 
equations  will  be  [41] 


ei,j  ~  Gj-i 


(23) 


Pi,j  =  Pi  J- 1  +  ar  (24) 

where  of  is  the  variance  of  the  process  noise  r.  Pi,j-i  is 
initialized  to  zero  for  the  first  sample  j  =  1.  The  measurement 
update  equations  will  be  [41] 


^ ,3 


Pi,  3 


kj  + 


h3 


(i_  _Ai. 


'h3 


2  kij  -  Aj) 

s 

)pi,j 


(25) 

(26) 


where  of  is  the  variance  of  the  measurement  noise  s.  Note 
that  since  both  the  state  and  measurement  z^j  are  scalars, 
all  the  update  equations  only  require  scalar  operations.  We 
estimate  of  and  of  as  following.  Let  the  mean  (calculated 
using  a  Gaussian  kernel  of  standard  deviation  of  gq  =  0.6 
sample)  of  the  strain  values  in  3  x  3  blocks  around  samples 
(i,j  —  1)  and  (i,j)  be  Pj-i  and  pj  respectively.  Then  of  is 
[52] 


ar  =  (P’j-l  ~  Pif-  (27) 

This  is  a  reasonable  estimate  of  of  as  it  tries  to  capture  the 
difference  between  pixel  values  at  adjacent  RF-lines.  If  the 
difference  between  the  mean  strain  values  is  high,  less  weight 
is  given  to  the  a  priori  estimate.  This  space-variant  estimation 
of  the  model  noise  provides  a  better  match  to  local  variations 
in  the  underlying  tissue  leading  to  a  greater  noise  reduction, 
of  is  the  variance  of  z^j  measurements  in  the  entire  image 
and  is  constant  throughout  the  image. 

The  strain  estimation  algorithm  can  be  summarized  as 
following: 

1)  Perform  least  squares  regression  in  the  axial  direction  for 
each  RF-line.  Generate  a  (noisy)  strain  image  Z  whose 
pixel  i,  j  is  z^j.  This  step  ensures  continuity  in  the  axial 
direction. 
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(a)  Unbiased  reg.  Entire  image 


depth  (mm) 

(b)  Schematic  displacements 


depth  (mm) 


(c)  Calculated  displacements  at  2%  strain 


(d)  Calculated  displacements  at  6%  strain  (e)  Unbiased  reg.  Top  of  the  image 


(f)  Biased  reg.  Entire  image 


Fig.  3.  Axial  strain  estimation  in  the  first  simulated  phantom,  (a)  The  SNR  values  corresponding  to  the  unbiased  regularization  calculated  in  the  entire  image, 
(b)  Schematic  plot  showing  the  underestimation  of  the  displacement  (Data  +  reg.  curve)  with  unbiased  regularization  (refer  to  the  text),  (c)-(d)  The  calculated 
displacements  at  the  bottom  of  a  RF-line  at  2%  strain  and  6%  strain  levels  respectively  with  biased  and  unbiased  regularization  terms.  The  ground  truth 
matches  the  displacement  given  by  the  biased  regularization  almost  perfectly,  and  therefore  is  not  shown  in  (c)  and  (d)  not  to  block  the  biased  regularization 
results.  The  length  of  the  RF-line  is  2560  (49.3  mm),  (e)  The  SNR  values  corresponding  to  the  unbiased  regularization  calculated  by  omitting  the  bottom  300 
samples  of  the  image,  (f)  The  SNR  values  corresponding  to  the  biased  regularization  calculated  in  the  entire  image.  Note  that  the  scale  of  the  SNR  in  graph 
(a)  is  much  smaller  than  that  of  graphs  (e)  and  (f). 


2)  Apply  the  Kalman  filter  to  the  noisy  strain  image  Z  in 
the  lateral  direction.  Generate  a  (denoised)  strain  image 
whose  pixel  i,  j  is  This  step  ensures  continuity  in 
the  lateral  direction. 


AM  method  according  to 


SNR  =  - 
a 


(28) 


Both  steps  are  applied  once  and  are  not  iterated.  We  show  in 
the  experimental  results  how  the  Kalman  filter  removes  the 
noise  from  the  strain  image  with  minimal  blurring,  owing  to 
the  model  noise  update  Equation  27. 


III.  Simulation  Results 

Field  II  [53]  and  ABAQUS  (Providence,  RI)  software  are 
used  for  ultrasound  simulation  and  for  finite  element  simu¬ 
lation.  Many  scatterers  are  distributed  in  a  volume  and  an 
ultrasound  image  is  created  by  convolving  all  scatterers  with 
the  point  spread  function  of  the  ultrasound  and  adding  the 
results  using  superposition.  The  phantom  is  then  meshed  and 
compressed  using  finite  element  simulation,  giving  the  3D 
displacement  of  each  node  of  the  mesh.  The  displacement 
of  each  scatterer  is  then  calculated  by  interpolating  the  dis¬ 
placement  of  its  neighboring  nodes.  Scatterers  are  then  moved 
accordingly  and  the  second  ultrasound  image  is  generated.  The 
displacement  and  strain  fields  are  then  calculated  using  the  AM 
methods  and  are  compared  with  the  ground  truth.  The  unitless 
metric  signal  to  noise  ratio  (SNR)  and  contrast  to  noise  ratio 
(CNR)  are  also  calculated  to  assess  the  performance  of  the 


where  st  and  are  the  spatial  strain  average  of  the  target  and 
background,  of  and  of  are  the  spatial  strain  variance  of  the 
target  and  background,  and  s  and  a  are  the  spatial  average 
and  variance  of  a  window  in  the  strain  image  respectively. 

The  parameters  of  the  ultrasound  probe  are  set  to  mimic 
commercial  probes.  The  probe  frequency  is  7.27  MHz,  the 
sampling  rate  is  40  MHz  and  the  fractional  bandwidth  is  60%. 
A  Hanning  window  is  used  for  apodization,  the  single  transmit 
focus  is  at  22.5  mm,  equi-distance  receive  foci  are  from  5  mm 
to  45  mm  at  each  5  mm,  the  transmit  is  sequential,  and  the 
number  of  active  elements  is  64. 

Two  simulated  phantoms  are  generated.  The  first  phantom 
is  50  x  10  x  55  mm  and  the  second  one  is  36  x  10  x  25mm. 
Respectively  5  x  105  and  1.4  x  105  scatterers  with  Gaussian 
scattering  strengths  [54]  are  uniformly  distributed  in  the  first 
and  second  phantom,  ensuring  more  than  10  scatterers  [55] 
exist  in  a  resolution  cell. 

The  mechanical  properties  of  both  phantoms,  required  for 
finite  element  simulation,  is  assumed  to  be  isotropic  and 
homogeneous.  The  first  phantom  is  uniform  while  the  second 
phantom  contains  a  circular  hole  filled  with  blood  that  can 
move  out-of-plane,  simulating  a  blood  vessel  in  tissue  (Figure 
7  (a)).  The  scatterers  are  distributed  in  the  vessel,  also  with  the 
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Fig.  4.  Lateral  strain  estimation  using  the  2D  AM  method  in  the  first 
simulated  phantom. 


same  intensity  and  distribution  as  the  surrounding  material.  A 
uniform  compression  in  the  z  direction  is  applied  and  the  3D 
displacement  field  of  phantoms  is  calculated  using  ABAQUS. 
The  Poisson’s  ratio  is  set  to  v  =  0.49  in  both  phantoms  to 
mimic  real  tissue  [56],  [57],  which  causes  the  phantoms  to 
deform  in  x  &  y  directions  as  a  result  of  the  compression  in 
the  z  direction. 

The  first  phantom  undergoes  compressions  in  the  z  direction 
to  achieve  strain  levels  of  1%  to  10%.  Figure  3  shows  the  SNR 
of  the  axial  strain  of  the  ID  AM  and  2D  AM  methods  (the 
window  for  SNR  calculation  covers  the  entire  strain  image  in 
(a)  &  (f)).  The  sharp  drop  of  the  SNR  with  strain  in  graph 
(a)  is  mainly  due  to  the  strain  underestimation  in  the  bottom 
part  of  the  image.  It  can  be  explained  as  following.  The  un¬ 
biased  regularization  term  tries  to  force  constant  displacement 
(dashed-dotted  red  line  in  (b)).  Assuming  an  ideal  noiseless 
case  where  the  data  term  gives  a  smooth  ramp  displacement 
(dashed  black  line  in  (b)),  minimizing  the  cost  function  (which 
is  the  summation  of  the  data  and  the  regularization  terms)  will 
underestimate  the  displacement  at  the  two  ends  (solid  blue 
line  in  (b)).  This  underestimation  decays  exponentially  moving 
towards  the  center  of  the  image.  This  artifact  is  shown  in  the 
simulation  experiment  at  2%  and  6%  strain  levels  in  (c)  and 
(d).  Since  we  exploit  the  fact  that  the  axial  displacement  of  the 
first  sample  is  zero  (Section  II-C),  the  underestimation  does 
not  happen  in  the  top  of  the  image.  Biasing  the  regularization 
prevents  this  artifact,  as  is  shown  in  (c)  and  (d).  The  AM 
method  with  or  without  the  bias  term  gives  the  same  result 
away  from  the  bottom  of  the  image:  part  (e)  shows  that  if  we 
ignore  300  (5.8  mm)  samples  at  the  bottom  of  the  image,  the 
SNR  will  not  drop  sharply  unlike  in  part  (a).  Part  (f)  shows  the 
SNR  of  the  AM  methods  with  biased  regularization  calculated 
in  the  entire  image.  The  SNR  at  1%  strain  in  parts  (e)  and  (f) 
is  the  same.  At  higher  strain  levels,  the  strain  underestimation 
propagates  more  into  the  middle  of  the  image,  and  therefore 
the  SNR  decreases  at  higher  strain  levels  in  graph  (e).  Part  (e) 
shows  2D  AM  gives  slightly  better  axial  strain  compared  to 
ID  AM.  IRLS  slightly  increases  the  SNR.  However,  we  will 
see  in  the  simulation  results  of  the  second  phantom  that  in 
the  presence  of  outliers  significant  improvement  in  SNR  and 
CNR  is  achieved  using  IRLS. 

The  SNR  of  the  lateral  strain  field  is  much  lower  than  that  of 
the  axial  strain  field  (Figure  4).  Unbiased  regularization  gives 
the  lowest  SNR,  mainly  due  to  artifacts  in  the  bottom  of  the 
image.  Similar  to  the  axial  strain,  the  SNR  improves  as  300 


(a)  Bias  (b)  Variance 

Fig.  5.  Bias  and  Variance  of  the  axial  strain  as  a  function  of  the  axial 
regularization  weight  a.  The  ground  truth  axial  and  lateral  strain  fields  are 
respectively  uniform  2%  and  2u%  fields  (y  =  0.49  is  the  Poisson’s  ratio).  The 
solid  blue  and  dashed  black  curves  both  correspond  to  unbiased  regularization 
and  the  solid  black  curve  corresponds  to  the  biased  regularization.  In  the  solid 
blue  and  solid  black  curves,  the  entire  image  is  included  in  the  calculation  of 
the  bias  and  noise.  In  the  dashed  black  curve  the  bottom  part  of  the  strain  field 
which  suffers  from  high  bias  (Figure  3  (b))  is  excluded  from  the  calculation  of 
the  bias  and  noise.  ID  AM  and  2D  AM  have  very  similar  bias  and  variance. 
The  curves  with  and  without  IRLS  are  also  very  close.  Therefore  each  curve 
corresponds  to  ID  AM  or  2D  AM  with  or  without  IRLS. 


samples  from  the  bottom  of  image  are  omitted  from  the  SNR 
calculation  (results  not  shown). 

The  effect  of  the  regularization  weights  on  bias  and  variance 
of  the  axial  strain  image  at  2%  ground  truth  axial  strain 
is  shown  in  Figure  5.  The  blue  curves  show  the  bias  and 
variance  of  the  entire  strain  image  obtained  with  unbiased 
regularization.  It  shows  the  tradeoff  between  the  bias  and 
variance:  increasing  the  regularization  weight  increases  the 
bias  and  decreases  the  variance.  The  variance  starts  to  increase 
at  a  ~  12  which  is  caused  by  the  underestimation  of  the  strain 
at  the  bottom  of  the  image  (the  artifact  in  Figure  3  (c)).  If  we 
exclude  the  bottom  300  samples  of  the  strain  image  from  the 
bias  and  variance  calculation  (the  black  dashed  curve),  we  see 
a  consistent  drop  of  variance  as  a  is  increased.  The  black 
curves  show  the  bias  and  variance  of  the  entire  strain  image 
obtained  with  biased  regularization.  Biasing  the  regularization 
causes  the  bias  to  decrease  as  the  regularization  weight  a  is 
increased  which  is  a  nonstandard  behavior.  It  can  be  explained 
by  the  simple  ground  truth  strain  field  which  is  uniform, 
exactly  what  the  regularization  term  is  trying  to  achieve.  Even 
in  the  unbiased  case,  only  the  bias  of  the  bottom  part  of  the 
strain  field  increases  as  a  is  increased  (i.e.  in  the  bias  plot,  the 
blue  curve  increases  while  the  black  dashed  curve  decreases). 
Therefore,  one  cannot  conclude  from  this  experiment  that 
higher  a  is  beneficial  to  both  bias  and  variance.  To  prove  this, 
we  designed  a  simulation  study  where  the  underlying  axial 
strain  field  continuously  varied  with  depth  and  the  lateral  and 
elevational  strains  were  zero  (such  strain  field  is  not  physically 
realizable).  We  observed  that  the  absolute  value  of  the  bias 
monotonically  increases  with  a  with  both  unbiased  and  biased 
regularizations.  To  save  space,  we  do  not  present  the  full 
results  here.  Similar  curves  for  the  lateral  strain  field  is  shown 
in  Figure  6. 

The  second  simulation  experiment  is  designed  to  show  the 
effect  of  smoothness  weight  and  IRLS  threshold  CNR  when 
the  correlation  is  lower  in  parts  of  the  image  due  to  fluid 
motion.  The  phantom  contains  a  vein  oriented  perpendicular 
to  the  image  plane  (Figure  7).  The  background  window  for 
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Fig.  7.  Second  simulation  experiment.  Measurements  in  (a)  are  in  mm.  In  (b),  a  scatterer  is  shown  in  the  bottom  left  part  as  a  red  dot.  Its  displacement  is 
calculated  by  interpolating  the  displacements  of  its  3  neighboring  nodes  on  the  mesh.  The  target  (circular)  and  background  (rectangular)  windows  for  CNR 
calculation  of  (d)  are  shown  in  (c). 


(a)  Bias  (b)  Variance 

Fig.  6.  Bias  and  Variance  of  the  lateral  strain  as  a  function  of  the  axial 
regularization  weight  a.  The  ground  truth  axial  and  lateral  strain  fields  are 
respectively  uniform  2%  and  2u%  fields  ( v  =  0.49  is  the  Poisson’s  ratio). 
The  solid  blue  curve  corresponds  to  unbiased  regularization  and  the  dashed 
and  solid  black  curves  correspond  to  the  biased  regularization.  IRLS  is  not 
used  in  the  solid  blue  and  dashed  black  curves. 


CNR  calculation  is  located  close  to  the  target  window  to  show 
how  fast  the  strain  is  allowed  to  vary,  a  property  related  to 
the  spatial  resolution.  The  maximum  CNR  with  IRLS  is  5.3 
generated  at  T  =  0.005  and  aa  =  38,  and  without  IRLS  is 
4.8  at  a a  =  338.  Such  high  aa  value  makes  the  share  of  the 
data  term  in  the  cost  function  very  small  and  causes  over¬ 
smoothing. 

A.  Displacement  Simulation 

To  study  the  performance  of  the  Kalman  filter,  we  simulate 
a  displacement  field  of  size  100  x  100  samples  whose  strain 
image  (calculated  using  least  squares  regression)  is  as  shown 
in  Figure  8  (a).  100  samples  in  the  axial  direction  corresponds 
to  approximately  1.9  mm  (assuming  40  MHz  sampling  rate), 
and  100  samples  in  the  lateral  direction  corresponds  to  10 
mm  to  25  mm  depending  on  the  probe.  To  be  consistent  with 
the  notations  of  Section  II-D,  let  denote  the  strain  values 
of  the  uncontaminated  image  in  (a).  We  then  contaminate 
the  displacement  field  with  a  Gaussian  noise  with  standard 
deviation  of  1.5  samples,  and  perform  least  squares  regression 
to  calculate  the  noisy  estimates  Zij  (Figure  8  (b)).  We  then 
apply  the  Kalman  filter  as  described  in  Section  II-D  to  the 


noisy  estimates  z^j  in  the  lateral  direction  (i.e.  row-by-row). 
The  posterior  estimates  of  the  strain  values,  are  shown  in 
(c).  The  strain  values  of  the  shown  line  in  (a)  -  (c)  (at  i  =  50 
samples)  is  shown  in  (d)  and  (e)  (The  plot  in  (d)  around  the 
step  in  magnified  in  (e)).  The  Kalman  filter  formulation  is 
eliminating  the  noise  without  over-smoothing  the  strain  image. 
This  is  due  to  the  model  variance  update  Equation  27.  We  note 
that  although  displacement  is  generally  continuous  in  tissue,  its 
spatial  derivation  (strain)  is  not:  at  the  boundary  of  two  tissues 
with  different  elasticity  moduli,  strain  field  is  discontinuous. 

IV.  Experimental  Results 

For  experimental  evaluation,  RF  data  is  acquired  from 
an  Antares  Siemens  system  (Issaquah,  WA)  at  the  center 
frequency  of  6.67  MHz  with  a  VF10-5  linear  array  at  a 
sampling  rate  of  40  MHz.  Only  the  2D  AM  method  is  used 
in  the  experimental  results.  Phantom  results  and  patient  trials 
are  presented  in  this  section.  The  tunable  parameters  of  the 
2D  AM  algorithm  are  set  to  a  =  5,  /3a  =  10,  Pi  =  0.005  and 
T  =  0.2  (Equations  12  &  20),  and  the  tunable  parameters  of 
the  DP  (run  for  the  seed  RF-line  in  the  2D  AM  algorithm) 
are  aa  =  ai  =  0.15  (Equation  1)  in  all  the  phantom  results 
(except  if  specified  otherwise).  In  the  patient  results,  all  the 
parameters  are  the  same  except  for  (3a  which  is  increased  to 
Pa  =  20  because  the  data  is  noisier.  The  strain  images  in  all 
the  patient  trials  are  obtained  using  the  least  squares  regression 
and  Kalman  filtering  as  described  in  Section  II-D. 

A.  Phantom  Results 

1 )  Effect  of  Regularization  on  Residuals:  The  cost  function 
of  the  AM  method  (Equation  7)  is  composed  of  residuals  (i.e. 
the  data  term)  and  the  regularization  terms.  The  AM  method 
minimizes  this  summation.  Therefore  the  AM  method  will  not 
necessarily  minimize  the  residuals.  We  now  show  that  the 
data  term  alone  is  non-convex  and  has  many  local  minima. 
Adding  the  regularization  term  will  eliminate  many  of  the  local 
minima  and  makes  optimization  of  the  data  term  easier.  This 
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Fig.  8.  Performance  of  the  least  squares  regression  and  the  Kalman  filter  on  simulated  motion  data,  (a)  shows  the  strain  field  calculated  using  least  squares 
regression  of  the  uncontaminated  displacement  field,  (b)  depicts  the  strain  field  calculated  using  least  squares  regression  of  the  contaminated  displacement 
field,  (c)  shows  the  strain  field  calculated  from  the  noisy  measurements  of  (b)  using  the  proposed  Kalman  filter  (KF  in  (b)  and  (c)  refers  to  Kalman  filter). 
The  pixels  of  images  in  (a)  to  (c)  are  respectively  the  ground  truth  (unavailable)  strain  values  eij,  the  noisy  measurements  Zij,  and  posterior  strain  values 
iij.  The  brightness  scale  in  (a)  to  (c)  is  the  same,  (d)-(e)  are  the  strain  estimation  at  the  horizontal  line  shown  in  (a)  to  (c).  (d)  is  magnified  in  (e)  around 
the  step.  The  Kalman  filter  removes  the  noise  while  keeping  the  image  sharp,  due  to  the  variable  model  noise  of  Equation  27. 


is  in  addition  to  the  effect  of  regularization  that  makes  the 
displacement  field  smooth,  a  generally  desired  attribute. 

The  effect  of  regularization  on  the  residuals  is  studied 
using  experimental  data.  An  elastography  phantom  (CIRS 
elastography  phantom,  Norfolk,  VA)  is  compressed  0.2  in 
axially  using  a  linear  stage,  resulting  in  an  average  strain  of 
6%.  Two  RF  frames  are  acquired  corresponding  to  before  and 
after  the  compression.  The  Young’s  elasticity  modulus  of  the 
background  and  the  lesion  under  compression  are  respectively 
33  kPa  and  56  kPa.  The  displacement  map  is  calculated  using 
the  2D  AM  method  and  the  residuals  corresponding  to  the 
displacement  map  are  obtained.  Figure  9  (a)-(c)  shows  the 
axial  and  lateral  strains  at  such  a  high  strain  rate  (minimum 
of  2%  and  maximum  of  11%).  The  mean  and  median  of 
the  residuals  p(ri)  in  the  entire  image  is  shown  in  (d). 
One  could  expect  the  graph  to  monotonically  increase  as 
the  regularization  weight  a  increases,  since  the  difference 
between  the  objective  function  C  and  the  residuals  E^1p(r,i) 
is  increased  as  a  is  increased.  However,  the  residual  values  are 
very  high  at  very  low  a.  Therefore,  numerical  minimization 
of  TifL1p{ri)  +  i2(Ad)  gives  a  smaller  value  for  Y^LiPiTi) 
compared  to  trying  to  directly  minimize  E^1p(r^).  This 
indicates  that  the  nonregularized  cost  function  is  not  quasi- 
convex  and  is  very  hard  to  minimize. 

2)  Resolution  of  the  Strain  Images  Generated  with  AM:  The 
effect  of  the  regularization  on  spatial  resolution  is  evaluated 
experimentally  using  the  experimental  setup  of  the  previous 
experiment.  The  compression  is  set  to  0.1  in  in  this  experiment. 


Figure  10  (a)  shows  the  strain  image  obtained  by  compression 
the  lesion  with  the  Young’s  modulus  of  56  kPa.  Spatial  reso¬ 
lution  is  evaluated  using  modulation  transfer  function  (MTF), 
an  established  method  for  estimating  the  spatial  resolution  of 
medical  imaging  systems  that  was  relatively  recently  extended 
to  elastography  [58].  The  spatial  resolution  of  the  recon¬ 
structed  images  is  determined  with  a  three-step  approach  [59], 
[60];  first,  the  edge  spread  function  is  computed  by  averaging 
the  pixel  values  across  the  background-inclusion  interface  (the 
line  in  Figure  10  (a));  second,  the  line  spread  function  (LSF) 
is  computed  by  differentiating  the  edge  spread  function;  third, 
the  MTF  is  determined  by  computing  the  Fourier  transform  of 
the  LSF  and  normalizing  the  resulting  function  to  zero  spatial 
frequency: 

MTF(fc)  =  (29) 

Figure  10  (c)  shows  the  MTF  for  five  different  normalization 
coefficients  respectively.  Strain  results  are  obtained  with  a 
regression  window  of  length  2k  +  1  =  65  (Section  II-D). 
Increasing  the  regularization  weight  is  adversely  affecting 
spatial  resolution.  Spatial  resolution  is  defined  as  the  spatial 
frequency  when  the  value  of  MTF  is  0.1.  At  a  =  1,  a  =  2 
and  a  =  4  this  value  is  respectively  2  cycles/mm,  1  cycles/mm 
and  0.5  cycles/mm.  In  addition  to  a,  this  value  also  depends 
on  the  length  of  the  regression  window  2k  +  1. 

3)  Image  Quality  Versus  Axial  and  Lateral  Sampling  Rates 
of  the  RF-Data:  Sampling  rate  of  the  RF-data  usually  ranges 
from  20  MHz  to  50  MHz  depending  on  the  hardware  of  the 
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Fig.  9.  Phantom  experimental  results.  The  top  row  shows  axial  displacement  and  axial  strains  as  labeled  (KF  in  (c)  refers  to  Kalman  filter).  Average  axial 
strain  and  maximum  strain  are  approximately  6.6%  and  11%.  (d)  and  (e)  show  lateral  displacement  and  lateral  strain  respectively,  (f)  shows  residuals  as  the 
regularization  weight  varies. 
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Fig.  10.  Phantom  experimental  results  showing  the  resolution  of  the  2D  AM.  (a)  Strain  image.  The  edge  spread  function  is  evaluated  along  the  vertical  line, 
(b)  The  strain  across  the  edge  (vertical  line  in  (a))  for  the  five  shown  regularization  values,  (c)  The  MTF  calculated  across  the  vertical  line  in  (a).  Spatial 
resolution  is  defined  as  the  spatial  frequency  when  the  value  of  MTF  is  0.1. 


device.  The  number  of  the  A-lines  provided  in  an  image  also 
varies  significantly.  In  addition,  bandwidth  limitations  of  the 
data  transfer  can  impose  limits  on  the  size  of  the  image  for 
real-time  operations.  In  this  study,  we  downsample  the  RF- 
data  by  a  factor  of  2  to  4  in  the  axial  direction  and  by  a 
factor  of  2  to  8  in  the  lateral  direction.  Figure  11  (a)-(g) 
shows  axial  and  lateral  displacement  and  strain  images  of  the 
CIRS  elastography  phantom  undergoing  maximum  axial  strain 
of  5%.  Axial  sampling  rate  can  be  reduced  by  a  factor  of  2 
without  significant  impact  on  the  strain  image  quality  (part 
(h)).  Downsampling  the  images  in  the  lateral  direction  by  a 
factor  of  4  results  the  CNR  of  the  axial  and  lateral  strain 
images  to  drop  respectively  12%  (from  16.3  to  14.3)  and  56% 
(from  2.55  to  1.13)  as  shown  in  (i).  While  the  axial  strain  is 
robust  to  the  number  of  A-line  in  the  image  even  at  a  high 


strain  level  of  5%,  the  lateral  strain  is  sensitive  to  it  (i).  Similar 
study  with  lower  axial  strain  levels  shows  that  as  the  axial 
strain  decreases,  higher  downsampling  rates  in  both  axial  and 
lateral  directions  are  possible  without  a  large  impact  on  the 
results. 

4)  Kalman  Filter:  The  performance  of  the  Kalman  filter 
is  studied  using  the  RF-data  used  in  Figure  9.  The  linear 
least  squares  differentiation  technique  is  applied  to  the  axial 
displacement  field  calculated  with  2D  AM,  resulting  in  Zij 
(Figure  12  (a)).  The  Kalman  filter  is  then  applied  to  Zij 
measurements  of  (a),  giving  the  posterior  measurements 
of  (b).  Comparing  the  strain  values  at  a  horizontal  line  of 
(a)  and  (b),  the  noisy  Zij  measurements  are  smoothed  in 
the  lateral  direction  using  the  proposed  Kalman  filter,  with 
minimal  blurring  of  the  edge. 
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Fig.  11.  Results  of  the  CIRS  elastography  phantom  at  5%  maximum  strain  at  different  axial  and  lateral  sampling  rates.  The  hard  lesion  is  spherical  and  has 
a  diameter  of  1  cm.  Downsampling  is  performed  by  simply  skipping  samples  in  the  axial  or  (and)  lateral  directions.  In  (c)  &  (f),  a  downsampling  ratio  of  2 
is  applied  in  both  axial  and  lateral  directions.  The  lateral  displacement  is  shown  in  number  of  samples  in  (d)  to  (f).  (h)  and  (i)  show  the  CNR  between  the 
target  and  background  windows  in  the  strain  images  as  the  axial  or  lateral  downsampling  rates  change.  The  target  and  background  windows  are  shown  in  the 
axial  strain  images  (a)  to  (c)  and  the  lateral  strain  image  (g).  In  (i),  the  lateral  strain  curve  is  not  calculated  for  downsampling  ratios  of  6  and  higher  because 
the  background  window  moves  out  of  the  image.  The  black  dashed  curve  with  the  highest  CNR  is  the  strain  obtained  with  the  Kalman  filter  (KF). 


B.  Clinical  Study 

Seven  patients  undergoing  open  surgical  radiofrequency 
(RF)  thermal  ablation  for  primary  or  secondary  liver  cancer 
were  enrolled  between  February  06,  2008  and  July  28,  2009. 
All  patients  enrolled  in  the  study  had  unresectable  disease 
and  were  candidates  for  RF  ablation  following  review  at 
our  institutional  multidisciplinary  conference.  Patients  with 
cirrhosis  or  suboptimal  tumor  location  were  excluded  from 
the  study.  All  patients  provided  informed  consent  as  part  of 
the  protocol,  which  was  approved  by  the  institutional  review 
board.  RF  ablation  was  administered  using  the  RITA  Model 
1500  XRF  generator  (Rita  Medical  Systems,  Fremont,  CA). 
Strain  images  are  generated  offline.  Some  preliminary  results 
are  published  in  [15]. 

We  show  the  results  from  only  4  patients  due  to  space 
limitations.  Figure  13  shows  the  B-mode  scan,  the  strain 


images  and  CT  scans  performed  after  RF  ablation.  Tissue  is 
simply  compressed  freehand  at  a  frequency  of  approximately 
1  compression  per  2  sec  with  the  ultrasound  probe  without 
any  attachment.  The  shadow  in  Figure  13  (a)  at  20  mm  depth 
is  produced  by  the  thermal  lesion.  Note  that  it  is  not  possible 
to  ascertain  the  size  and  position  of  the  thermal  lesions  from 
B-mode  images.  In  addition,  the  thermal  lesion  has  different 
appearances  in  the  three  B -scans.  However,  the  thermal  lesions 
show  very  well  as  hard  lesions  in  the  strain  images.  After 
gross  correlation  of  the  post  ablation  CT  scan  and  the  thermal 
lesion  in  the  strain  images,  the  size  of  the  lesion  seems  to 
correspond  well.  However,  a  more  rigorous  validation  of  the 
size  and  shape  of  the  ablated  lesion  in  the  elastography  image 
is  underway  using  nonrigid  registration  of  CT  and  ultrasound 
images.  To  the  best  of  our  knowledge,  this  is  also  the  first 
demonstration  of  the  success  of  elastography  in  imaging  the 
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Fig.  12.  Performance  of  the  least  squares  regression  and  the  Kalman  filter  in  experimental  data,  (a)  shows  the  axial  strain  field  calculated  by  least  squares 
regression  of  the  noisy  displacement  field,  (b)  depicts  the  strain  field  calculated  from  the  noisy  measurements  of  (a)  using  the  proposed  Kalman  filter  (KF  in 
(a)  and  (b)  refers  to  Kalman  filter).  The  pixels  of  images  in  (a)  and  (b)  are  respectively  the  least  squares  measurements  Zij,  and  posterior  strain  values  eij. 
(c)  shows  the  strain  estimation  at  the  17  mm  deep  horizontal  line  shown  in  (a)  to  (b).  The  Kalman  filter  removes  the  noise  while  keeping  the  image  sharp, 
due  to  the  variable  model  noise  of  Equation  27. 


thermal  lesion  in  an  in-vivo  human  experiment. 

We  have  also  acquired  patient  RF  data  of  liver  ablation 
prior  and  after  ablation  in  one  of  the  patient  trials.  Figure 
14  shows  the  B-mode,  strain  and  venous  and  arterial  phase  2 
CT  images  obtained  before  ablation,  and  Figure  15  shows  the 
B-mode,  strain  and  lateral  displacement  images  after  ablation. 
In  Figure  14,  the  tumor  (marked  in  the  CT  images  (f)  and  (g)) 
is  not  visible  in  the  B-mode  image  (a),  but  is  clearly  visible 
in  the  strain  images  (b)  and  (c).  While  the  tissue  is  getting 
compressed  with  the  ultrasound  probe,  the  middle  hepatic  vein 
(marked  as  5)  which  is  only  4-8  cm  from  vena  cava  inferior 
pulsates  at  high  amplitude.  The  graph  in  (e)  schematically 
shows  the  probe  motion  and  variations  in  the  diameter  of  the 
vein.  Therefore,  the  vein  can  look  soft  as  in  (c)  or  hard  as  in 
(b)  depending  on  whether  its  diameter  variation  is  in  the  same 
(marked  by  ellipse  1  in  (e))  or  opposite  (marked  by  ellipse  2 
in  (e))  direction  as  the  probe  motion.  The  effect  of  pulsation 
of  vessels,  a  well-known  cause  of  signal  decorrelation,  is 
minimized  via  IRLS  resulting  in  a  low  noise  strain  image. 
In  addition,  since  the  2D  AM  method  gives  a  dense  motion 
field  (same  size  as  RF  data),  the  small  artery  at  the  diameter 
of  less  than  2  mm  (marked  as  4  in  (a))  is  discernible  in  (b) 
from  the  low  pressure  portal  vein.  The  ablated  lesion  is  also 
discernible  in  the  strain  images  of  Figure  15  (b)  and  (c).  We 
believe  the  soft  region  in  the  middle  of  the  two  hard  ablation 
lesions  in  (b)  and  (c)  (at  the  depth  of  25-30  mm  and  width  of 
10-25mm)  is  not  close  to  any  of  the  10  tines  of  the  ablation 
probe.  Therefore  because  of  its  proximity  to  veins  and  vessels 
its  temperature  has  remained  low. 

V.  Discussion 

The  resolution  of  the  method  is  formally  studied  in  Section 
IV-A  using  the  phantom  experiment.  Future  work  will  include 
more  intuitive  measures  for  resolution  in  terms  of  the  smallest 
detectable  target  as  a  function  of  its  elasticity  difference  with 
the  background. 

2CT  scans  are  performed  at  different  phases  after  intravenous  injection  of 
a  contrast  agent.  In  the  arterial  phase  (directly  after  injection  of  a  contrast 
agent)  arteries  will  enhance,  where  as  in  the  venous  phase  (30-60  sec  after 
injection)  the  hepatic  parenchyma  and  veins  will  enhance. 


The  cost  function  is  a  regularized  function  of  all  dis¬ 
placements  on  an  A-line.  This  makes  the  methods  robust  to 
noise  which  exist  throughout  the  image.  Besides,  the  AM 
methods  are  not  window-based  and  therefore  they  don’t  suffer 
from  decorrelation  within  the  window.  As  a  result,  both  AM 
methods  work  for  strains  as  high  as  10%.  In  addition,  the  IRLS 
outlier  rejection  technique  makes  the  AM  methods  robust  to 
local  sources  of  decorrelation  such  as  out-of-plane  motion  of 
movable  structures  or  blood  flow. 

Global  stretching  assumes  a  constant  strain  across  the  depth 
and  stretches  one  of  the  RF-limes  accordingly.  It  is  shown 
that  it  enhances  the  quality  of  correlation  based  elastography 
methods.  The  reason  is  that  the  strain  of  each  point  can  be 
assumed  to  be  the  global  strain  (fixed  for  each  RF-line)  plus 
some  perturbation,  i.e.  constant  strain  is  a  better  approximation 
than  zero  strain.  Biasing  the  regularization  is  motivated  by  the 
same  reason  and  involves  almost  no  additional  computational 
cost. 

Improvement  in  the  SNR  and  CNR  achieved  with  Kalman 
filtering  differentiation  is  due  to  utilizing  the  (piecewise) 
continuity  of  the  strain  field.  One  could  think  of  a  unified 
framework  which  includes  both  the  2D  AM  and  the  Kalman 
filtering  and  directly  calculates  the  strain  field.  We  made  an 
effort  to  formulate  Equation  15  in  terms  of  strain  values. 
Unfortunately,  the  coefficient  matrix  in  the  L.H.S.  became 
a  full  matrix  for  our  desired  regularization.  Such  large  full 
system  cannot  be  solved  in  real-time. 

The  least  squares  differentiation  of  Section  II-D  can  be 
incorporated  in  the  Kalman  filter.  This  can  be  simply  done 
by  defining  the  state  at  each  point  to  be  the  displacement 
and  the  strain  of  that  point.  The  observed  variables  are  the 
noisy  displacement  measurements  from  2D  AM.  Solving  for 
the  state  gives  a  strain  estimate  at  each  point.  However,  we 
preferred  to  follow  the  common  approach  of  first  finding  the 
strain  by  solving  least  squares.  In  addition,  the  axial  and 
lateral  displacements  can  be  considered  as  two  channels  of 
a  measurement  and  a  Kalman  filter  that  takes  into  account 
both  intra-channel  (spatial)  and  inter-channel  variations  can 
be  developed.  This  is  a  subject  of  future  work. 

Lateral  displacement  estimation  with  2D  AM  is  of  order  of 
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Fig.  13.  In-vivo  images  of  the  thermal  lesion  produced  by  RF  ablation  therapy  of  liver  cancer.  All  images  acquired  after  ablation.  1st,  2nd  and  3rd  rows 
correspond  to  the  1st,  2nd  and  3rd  patients  respectively.  The  thermal  lesion  shows  in  (b),  (f)  &  (j)  as  dark,  surrounded  by  normal  tissue  in  white.  The  lateral 
displacement  images  are  shown  in  number  of  samples  (they  do  not  immediately  carry  anatomical  information).  In  (b),  (d),  (f),  (h),  (j)  &(1)  the  delineated 
thermal  lesions  is  shown.  The  non-unity  aspect  ratio  in  the  axes  of  the  B-mode  and  strain  images  should  be  considered  when  comparing  them  to  the  CT  scans. 


magnitude  less  accurate  than  the  axial  displacement  estimates. 
We  tested  the  following  algorithm  for  calculating  the  lateral 
displacement  field  based  on  ID  AM:  run  ID  AM  to  find 
the  axial  displacement  field  A,  then  transpose  both  ultrasound 
images  I\  and  I 2  and  run  ID  AM  again  using  A  calculated  in 
the  previous  step.  The  axial  displacement  field  calculated  for 
the  transposed  images  is  in  fact  the  lateral  displacement  of  the 
original  images.  Although  considerably  more  computationally 
expensive  than  2D  AM,  this  algorithm  did  not  improve  the 
lateral  displacement  estimation.  Therefore  only  images  of 
lateral  displacement  are  provided  for  the  patient  trials  because 
the  lateral  strain  did  not  show  the  ablation  lesion.  This  is 
in  accordance  with  recent  work  [36]  which  only  shows  the 
lateral  displacement.  A  2D  displacement  field  can  be  utilized 
to  calculate  the  thermal  expansion  and  to  reconstruct  the  strain 
tensor.  Incorporation  of  the  synthetic  lateral  phase  [61],  [62], 
[63]  into  2D  AM  to  further  improve  the  accuracy  of  the  lateral 
displacement  measurement  is  also  a  subject  of  future  work. 

In  cases  where  the  two  ultrasound  frames  correlate  very 
poorly  throughout  the  image,  ID  AM  outperforms  2D  AM 
because  DP  is  run  for  the  entire  image  in  ID  AM.  However, 


in  those  cases  the  strain  images  are  of  very  low  quality  even 
with  ID  AM.  In  cases  where  the  images  correlate  reasonably, 
the  2D  AM  algorithm  slightly  outperforms  ID  AM  in  terms 
of  the  SNR  of  the  axial  strain  as  shown  in  Figure  3  (e)  and 
(f).  Also,  ID  AM  and  2D  AM  are  very  similar  in  terms  of 
bias  and  variance  as  mentioned  in  the  caption  of  the  Figure 
5.  And  finally,  2D  AM  is  more  than  10  times  faster  than 
ID  AM  because  it  eliminates  the  redundant  calculations  in 
the  DP  step  of  ID  AM.  This  is  important  considering  that 
there  are  combinatorial  many  ways  of  choosing  two  frames 
for  elastography  from  a  sequence  of  images.  Having  a  fast 
algorithm,  like  2D  AM,  makes  it  plausible  to  invest  time  to 
perform  real-time  frame  selection,  an  area  that  we  are  currently 
working  on  [16],  [64]. 

Recent  work  [65]  has  attempted  to  reconstruct  elasticity 
from  the  displacement  field  for  monitoring  thermal  ablation.  It 
has  also  shown  that  [66]  compared  to  strain  images,  elasticity 
images  have  both  higher  correlation  with  the  ablation  zone  and 
give  higher  CNR.  Another  work  [67]  has  utilized  the  solution 
of  the  elasticity  reconstruction  to  improve  motion  estimation  in 
an  iterative  framework.  Calculation  of  the  elasticity  modulus 
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Fig.  14.  In-vivo  images  of  the  fourth  patient  before  RF  ablation.  In  (a),  the  left  anterior  branch  of  portal  vein  is  marked  as  1  and  2  and  has  low  pressure  and 
therefore  compresses  easily.  Arteries  (marked  as  3  and  4)  and  the  middle  hepatic  vein  (marked  as  5)  however  pulsate  with  the  heart  beat  and  may  have  low 
or  high  pressure,  (b)  and  (c)  both  show  the  axial  strain  from  the  same  location  before  ablation.  They  are  calculated  at  two  different  phases  of  the  heart  beat. 
The  cancer  tumor  is  discernible  in  (b)  and  (c)  (regardless  of  the  systolic  or  diastolic  blood  pressure),  and  its  boundary  is  shown.  1  and  2  (as  marked  in  (a)) 
correspond  to  the  high  strain  area  in  both  (b)  and  (c).  Since  3,  4  and  5  (as  marked  in  (a))  pulsate,  they  may  look  hard  (as  in  (b))  or  soft  (as  in  (c)).  (d)  shows 
the  lateral  displacement.  The  tumors  is  not  visible  in  this  image,  (e)  shows  the  motion  of  the  probe  and  the  variation  in  the  diameter  of  the  arteries  due  to 
the  heart  beat  (refer  to  the  text),  (f)  is  the  arterial  phase  and  (g)  is  the  venous  phase  contrast  CT  images.  The  numbers  1-5  mark  the  same  anatomy  as  (a). 
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Fig.  15.  In-vivo  images  of  the  fourth  patient  after  RF  ablation.  Similar  to  Figure  14,  the  hepatic  vein  (marked  as  5)  can  have  low  strain  (as  in  (b))  or  high 
strain  (as  in  (c))  values. 


in  our  ablation  monitoring  trials  is  an  area  of  future  work. 

Statistical  analysis  of  the  residuals  is  a  subject  of  future 
work.  The  sum  of  squared  differences  used  as  the  similarity 
metric  in  our  cost  function  is  suitable  if  ultrasound  noise  can 
be  modeled  as  additive  Gaussian  noise.  However,  ultrasound 
noise  is  not  simply  additive  Gaussian  and  it  has  been  shown 
that  similarity  metrics  that  model  the  noise  process  consid¬ 
ering  physics  of  ultrasound  give  more  accurate  results  [68]. 
Performance  of  the  2D  AM  method  for  images  that  are  not 
fully  developed  speckles  (i.e  have  few  scatterers  per  resolution 
cell)  is  also  a  subject  of  future  work. 

Current  implementations  of  the  ID  AM  and  2D  AM  take 
respectively  0.4  sec  and  0.04  sec  to  generate  strain  images 
(axial  for  ID  AM  and  axial  and  lateral  for  2D  AM)  of  size 
1000  x  100  on  a  3.8  GHz  P4  CPU.  DP  contributes  to  more 
than  90%  of  the  running  time  of  the  ID  AM,  and  that’s  why 


it  is  slower  than  2D  AM  where  DP  is  only  run  for  a  single 
A-line.  The  running  time  of  both  methods  changes  linearly 
with  the  size  of  the  image. 

VI.  Conclusion 

Two  regularized  elastography  methods,  ID  AM  and  2D  AM, 
are  introduced  for  calculating  the  motion  field  between  two 
ultrasound  images.  They  both  give  dense  subsample  motion 
fields  (ID  AM  gives  subsample  axial  and  integer  sample  lateral 
and  2D  AM  gives  subsample  axial  and  lateral)  in  real-time. 
The  size  of  the  motion  fields  is  the  same  as  the  size  of  the 
RF-data  (except  for  few  samples  from  the  boundary  whose 
displacements  are  not  calculated).  Such  dense  motion  fields 
lead  to  dense  strain  fields  which  are  critical  in  detecting  small 
lesions.  The  prior  of  tissue  motion  continuity  is  exploited  in 
the  AM  methods  to  minimize  the  effect  of  signal  decorrelation. 
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The  regularization  term  is  biased  with  the  average  strain  in  the 
image  to  minimize  underestimation  of  the  strain  values.  Parts 
of  the  image  that  have  very  low  correlation  are  treated  as  out¬ 
liers  and  their  effect  is  minimized  via  IRLS.  The  strain  image 
is  calculated  by  differentiating  the  motion  fields  using  least 
squares  regression  and  Kalman  filtering.  The  performance  of 
the  proposed  elastography  algorithms  is  analyzed  using  Field 
II  and  finite  element  simulations,  and  phantom  experiments. 
Clinical  trials  of  monitoring  RF  ablation  therapy  for  liver 
cancer  in  four  patients  are  also  presented.  An  implementation 
of  the  2D  AM  method,  the  least  squares  regression  and  the 
Kalman  filter  in  MATLAB  mex  functions,  as  well  as  some 
of  the  phantom  and  patient  RF  data  used  in  this  work  are 
available  for  academic  research  and  can  be  downloaded  from 
http://www.cs.jhu.edu/~rivaz/Ultrasound_Elastography/. 
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Abstract.  Tracked  ultrasound  elastography  can  be  used  for  guidance  in 
partial  breast  radiotherapy  by  visualizing  the  hard  scar  tissue  around  the 
lumpectomy  cavity.  For  clinical  success,  the  elastography  method  needs 
to  be  robust  to  the  sources  of  decorrelation  between  ultrasound  images, 
specifically  fluid  motions  inside  the  cavity,  change  of  the  appearance  of 
speckles  caused  by  compression  or  physiologic  motions,  and  out-of-plane 
motion  of  the  probe.  In  this  paper,  we  present  a  novel  elastography  tech¬ 
nique  that  is  based  on  analytic  minimization  of  a  regularized  cost  func¬ 
tion.  The  cost  function  incorporates  similarity  of  RF  data  intensity  and 
displacement  continuity,  making  the  method  robust  to  small  decorre¬ 
lations  present  throughout  the  image.  We  also  exploit  techniques  from 
robust  statistics  to  make  the  method  resistant  to  large  decorrelations 
caused  by  sources  such  as  fluid  motion.  The  analytic  displacement  esti¬ 
mation  works  in  real-time.  Moreover,  the  tracked  data,  used  for  targeting 
the  radiotherapy,  is  exploited  for  discarding  frames  with  excessive  out- 
of-plane  motion.  Simulation,  phantom  and  patient  results  are  presented. 


1  Introduction 

Breast  irradiation  after  lumpectomy  significantly  reduces  the  risk  of  cancer  re¬ 
currence.  There  is  growing  evidence  suggesting  that  irradiation  of  only  the  in¬ 
volved  area  of  the  breast,  partial  breast  irradiation  (PBI),  is  as  effective  as  whole 
breast  irradiation  [1].  Benefits  of  PBI  include  significantly  shortened  treatment 
time  and  fewer  side  effects  as  less  tissue  is  treated.  However,  these  benefits  cannot 
be  realized  without  localization  of  the  lumpectomy  cavity.  Tracked  ultrasound 
elastography  can  be  used  for  localizing  the  lumpectomy  cavity  in  the  treatment 
room,  minimizing  tissue  motion  from  planning  to  treatment. 

This  paper  is  focused  on  freehand  palpation  elastography,  which  involves  es¬ 
timating  the  displacement  field  of  the  tissue  undergoing  slow  compression.  Most 
elastography  techniques  estimate  the  displacement  field  using  local  cross  corre¬ 
lation  analysis  of  echoes  [2,3,4].  These  methods  are  very  sensitive  and  accurate 
for  calculating  small  displacements.  However,  elastography  is  subject  to  speckle 
decorrelation  caused  by  various  sources  such  as  motion  of  subresolution  scatter¬ 
ed,  out-of-plane  motion,  high  compression  and  complex  fluid  motions. 

The  prior  of  tissue  deformation  continuity  can  be  used  to  make  elastography 
more  robust  to  signal  decorrelation.  Previous  work  on  regularized  elastography 
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is  computationally  expensive  [5,6].  Dynamic  programming  (DP)  can  be  used 
to  speed  the  optimization  procedure  [7],  but  it  only  gives  integer  displacements. 
Subpixel  displacement  estimation  is  possible  [7],  but  it  is  computationally  expen¬ 
sive  if  a  fine  subpixel  level  is  desired.  In  addition,  a  fixed  regularization  weight  is 
applied  throughout  the  image.  However,  while  two  ultrasound  images  may  cor¬ 
relate  well  in  most  parts,  they  can  have  small  correlation  in  specific  parts.  Four 
examples  of  low  correlation  are:  (1)  correlation  decreases  with  depth  mainly  due 
to  a  decrease  in  the  ultrasonic  signal  to  noise  ratio,  (2)  correlation  is  low  close 
to  arteries  due  to  complex  motion  and  inside  vessels  due  to  blood  motion,  (3) 
correlation  is  extremely  low  in  lesions  that  contain  liquid  due  to  the  incoherent 
fluid  motion  [8,3],  and  (4)  out-of-plane  motion  of  movable  structures  within  the 
image  [8]  causes  low  local  correlation.  To  prevent  such  regions  from  introduc¬ 
ing  errors  in  the  displacement  estimation  one  should  use  large  weights  for  the 
regularization  term,  resulting  in  over-smoothing. 

Freehand  palpation  elastography  provides  ease-of-use  and  requires  minimum 
additional  cost.  However,  out-of-plane  motion  cannot  be  avoided  in  freehand  pal¬ 
pation,  which  reduces  the  quality  of  any  elastography  method.  Assisted  freehand 
elastography  [9]  significantly  reduces  the  out-of-plane  motion  but  it  requires  ad¬ 
dition  of  a  device  to  the  probe.  Quality  metrics  such  as  persistence  in  strain 
images  have  also  been  developed  to  address  this  problem  [10].  To  measure  the 
persistence,  elastography  is  performed  on  two  pairs  of  images  and  the  resulting 
strain  images  are  correlated.  This  method  requires  strain  images  for  calculating 
the  quality  metric.  Therefore,  trying  all  the  combinations  in  a  series  of  frames 
to  find  the  best  pair  for  elastography  will  be  computationally  expensive. 

In  this  paper,  we  present  a  novel  elastography  method  based  on  analytic 
minimization  (AM)  of  a  cost  function  that  incorporates  similarity  of  echo  ampli¬ 
tudes  and  displacement  continuity.  We  introduce  a  novel  regularization  term  and 
demonstrate  that  it  minimizes  displacement  underestimation  caused  by  smooth¬ 
ness  constraint.  We  also  introduce  the  use  of  robust  statistics  implemented  via 
iterated  reweighted  least  squares  (IRLS)  to  treat  uncorrelated  ultrasound  data 
as  outliers.  And  finally,  we  use  the  tracking  information  to  select  the  best  pairs 
of  frames  for  elastography.  Simulation,  phantom  and  patient  experiments  are 
presented  for  validation. 

2  Regularized  Displacement  Estimation 

Dynamic  Programming  (DP).  DP  is  a  discrete  efficient  optimization  tech¬ 
nique  for  causal  systems.  In  DP  elastography  [7],  a  cost  function  is  defined  as 

C(i,di)  =  min  {C(i  -  l,d*-i)  +  aaR(di,  di-i)}  +  \h(i)  -/2(i  +  d*)| ,  i  =  2---ra 

di- 1 

(1) 

where  di  is  the  displacement  of  sample  i,  R{di,di-\)  =  (di  —  di- 1)  is  an  axial 
regularization  term  (axial,  lateral  and  out-of-plane  directions  are  respectively 
z,  x  and  y  in  Figure  2  (a)),  aa  is  a  weight  for  the  regularization,  I\  and  I2 
are  corresponding  RF-lines  of  before  and  after  deformation  and  m  is  the  length 
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of  RF-lines.  The  cost  function  is  minimized  at  i  =  m  and  the  di  values  that 
have  minimized  the  cost  function  are  traced  back  to  i  =  1,  giving  the  di  for  all 
samples.  We  have  implemented  a  2D  DP  algorithm  similar  to  [7]  to  generate 
integer  displacements  as  a  starting  point  for  the  next  step  of  our  algorithm. 

Analytic  Minimization  (AM).  We  now  propose  a  method  that  analytically 
minimizes  a  regularized  cost  function  and  gives  the  refined  displacement  field. 
Only  axial  displacements  will  be  refined  for  strain  calculation. 

Having  the  integer  displacements  di  from  DP,  it  is  desired  to  find  Adi  values 
such  that  di  +  Adi  gives  the  value  of  the  displacement  at  the  sample  i  for  i  = 
1  •  •  •  m.  Such  Adi  values  will  minimize  the  following  regularized  cost  function 

C  (Adi,  •  •  • ,  Adm )  =  r™!  [h(i)  -I2(i  +  di  +  Adi)}2  + 
aa(di  +  Adi  ~  di- $  -  Adi- 1)2  +  a/(d»  +  Adt  -  df  -  Ad\')2  (2) 

where  superscript  p.  refers  to  the  previous  RF-line  (adjacent  RF-line  in  the 
lateral  direction)  and  ai  is  a  weight  for  lateral  regularization.  Substituting  l2(i  + 
di~\~  Adi)  with  its  first  order  Taylor  expansion  approximation  around  di,  we  have 

C  {Adi,  •  •  • ,  Adm)  =  £™1  lh(i)  -  I2(i  +  di)  -  /'(*  +  di)Adi )f  + 
aa(di  +  Adi  ~  di-i  -  Adi- i)2  +  «/(<**  +  Adt  ~  df  -  Adf)2  (3) 

where  J2  is  the  derivative  of  the  I2 .  The  optimal  Adi  values  occur  when  the 
partial  derivative  of  C  w.r.t.  Adi  is  zero.  Setting  dd^d  =  0  we  have 


(I'22  +  aaT>  +  ati)Ad  =  I'2e  -  (aaD  +  azi)d  +  ,  D  = 


1  -1  0  •••  0 

-1  2  -1  •••  0 


[0  •••  0  -1  1J 

(4) 

where  I2  =  diag(/2(l  +  di)  •  •  •  J2(ra  +  dm)),  Ad=  [Adi  •  •  •  Adm]  ,  e=  [e\  •  •  •  em]  , 
ei  =  I\(i)  —  12(1  +  di),  d  =  [di  *  •  •  dm]T,  d t-p-  =  dp  +  Adp  is  the  vector  of  total 
displacement  of  the  previous  line  and  I  is  the  identity  matrix.  I2,  D  and  I  are 
matrices  of  size  m  x  m  and  Ad ,  r,  d  and  dt  p-  are  vectors  of  size  m. 


Biasing  the  Regularization.  The  regularization  term  aa(di  +  Adi  —  di-i  ~ 
Adi- 1)2  penalizes  the  difference  between  di~\~Adi  and  di-\+Adi-\,  and  therefore 
can  result  in  underestimation  of  the  displacement  field.  Such  underestimation 
can  be  prevented  by  biasing  the  regularization  by  e  to  aa(di  +  Adi  ~  di-i  — 
Adi- 1  —  e)2,  where  e  =  (dm  —  d\)/{m  —  1)  is  the  average  displacement  difference 
between  samples  i  and  i  —  1.  An  accurate  enough  estimate  of  dm  —  d\  is  known 
from  the  previous  line.  With  the  bias  term,  the  R.H.S.  of  Equation  4  becomes 
I2e  —  (<aaD  +  ail)d  +  cqdt,;p-  +  b  where  the  bias  term  is  b  =  aa[— e  0  •  •  •  0  e]T 
and  all  other  terms  are  as  before.  Interestingly,  except  for  the  first  and  the  last 
equation  in  this  system,  all  other  m  —  2  equations  are  same  as  Equation  4. 

Equation  4  can  be  solved  for  Ad  in  4m  operations  since  the  coefficient  matrix 
I22  +  aaT>  -\-aiI  is  tridiagonal.  Utilizing  its  symmetry,  the  number  of  operations 
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can  be  reduced  to  2m.  The  number  of  operations  required  for  solving  a  system 
with  a  full  coefficient  matrix  is  more  than  m3/ 3,  significantly  more  than  2m. 

Making  Tracking  Resistant  to  Outliers.  Even  with  pure  axial  compression, 
some  regions  of  the  image  may  move  out  of  the  imaging  plane  and  increase 
the  decorrelation.  In  such  parts  the  confidence  of  the  data  term  is  less  and 
therefore  the  weight  of  the  regularization  term  should  be  increased.  The  parts 
of  the  image  with  low  correlation  can  be  regarded  as  outliers  and  therefore 
a  robust  estimation  technique  can  limit  their  effect.  Before  deriving  a  robust 
estimator  for  Z\d,  we  rewrite  Equation  3  as  C(Ad)  =  T’]T1p(r^)  +  R(Ad)  where 
ri  =Ii{i)  —  l2{iJrdi)  —  H2{iJtdi)Adi1p{ri)  =  rf  and  R  is  the  regularization  term. 
The  M-estimate  of  Ad  is  Ad  =  arg ruined  {T’]T1p(r^)  +  R(Ad)}  where  p(u)  is 
a  robust  loss  function  [11].  The  minimization  is  solved  by  setting  dd^d:  =  0: 


p\n) 


dr 

dAdi 


dR(Ad) 

dAdi 


(5) 


A  common  next  step  [11]  is  to  introduce  a  weight  function  re,  where  w(ri).ri  = 
p'(ri).  This  leads  to  a  process  known  as  “iteratively  reweighted  least  squares” 
(IRLS),  which  alternates  steps  of  calculating  weights  w(ri)  for  r%  =  1  •  •  •  m  using 
the  current  estimate  of  Ad  and  solving  Equation  5  to  estimate  a  new  Ad  with 
the  weights  fixed.  Among  many  proposed  shapes  for  w(-),  we  use  [11] 


w(ri) 


i  \n\<T 

w\  N>T 


(6) 


where  T  is  a  threshold  that  can  be  tuned.  A  small  T  will  treat  many  samples  as 
outliers.  With  the  addition  of  the  weight  function,  Equation  5  becomes 


(wl^2  +  aD  +  a2l)Ad  =  wl^e  —  (oqD  +  o^ijd  +  a2dt'p'  +  b  (7) 


where  w  =  diag(w(r\)  •  •  -w(rm)).  All  of  the  results  presented  in  this  work  are 
obtained  with  one  iteration  of  the  above  equation  unless  otherwise  specified. 
Current  implementation  of  the  AM  algorithm  with  the  IRLS  takes  0.015s  to 
generate  a  dense  displacement  field  of  size  1300  x  60  on  a  3.4GHz  P4  CPU(not 
including  the  DP  run  time).  The  computation  time  increases  linearly  with  the 
size  of  images. 

Frame  Selection.  The  ultrasound  probe  is  tracked  in  navigation/guidance 
systems  to  provide  spatial  information,  to  generate  freehand  3D  ultrasound, 
or  to  facilitate  multi-modality  registration.  Through  a  calibration  process,  the 
6DOF  motion  of  the  probe  in  the  sensor  coordinate  system  is  transformed  into 
image  coordinate  system  [12].  The  mean  of  the  absolute  motion  value  of  all  pixels 
in  3D,  (1^1),  {\vy\)  and  (\vz\),  can  be  analytically  related  to  the  6DOF  sensor 
readings  using  straightforward  and  efficient  geometric  computations.  For  frame 
i  and  j  to  be  selected  from  a  sequence  of  frames  for  elastography, 


Qi 


kx  (|^|)2  +  ky  ( \vy\ )2  +  kz 


<N> 


Uz,opt  | 


(N>  +  c 


(8) 
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should  be  minimized  where  kx,  ky,  and  kz  are  weights  for  lateral,  out-of-plane 
and  axial  displacements  and  vz^opt  is  the  optimum  axial  motion.  Please  refer  to 
[12]  for  a  rationale  of  the  shape  this  function.  Note  that  the  selected  pairs  are 
not  necessarily  consecutive  frames.  The  parameters,  kx,  ky,  kz ,  vZj0pt  and  c  are 
manually  tuned  to  1,  2,  1,  0.7  and  1  for  the  AM  elastography  method. 

3  Simulation,  Phantom  and  Patient  Results 

Simulation  Results.  RF  ultrasound  data  of  two  phantoms  are  simulated  using 
Field  II  [13].  The  first  phantom  is  50  x  10  x  55mm  and  the  second  one  is  36  x 
10  x  25mm.  They  are  both  made  of  homogeneous  and  isotropic  material:  the 
first  one  is  uniform  and  the  second  one  contains  a  circular  hole  filled  with  water, 
simulating  a  blood  vessel  in  tissue  (Figure  2  (a)).  A  uniform  compression  in  the 
z  direction  is  applied  and  the  3D  displacement  field  of  the  phantom  is  calculated 
using  ABAQUS  finite  element  package  (Providence,  RI).  The  Poisson’s  ratio  is 
set  to  v  =  0.49  in  both  phantoms  to  mimic  real  tissue,  which  causes  the  phantom 
to  deform  in  x  &  y  directions  as  a  result  of  the  compression  in  the  z  direction. 

Respectively  5  x  105  and  1.4  x  105  scatterers  with  uniform  scattering  strengths 
are  uniformly  distributed  in  the  first  and  second  phantom,  ensuring  more  than 
10  scatterers  exist  in  a  resolution  cell.  The  scatterers  are  distributed  in  the  8mm 
diameter  vein  also  (Figure  2  (a)).  To  construct  deformed  ultrasound  images,  the 
displacement  of  all  of  the  scatterers  is  calculated  by  interpolating  the  displace¬ 
ment  of  the  neighboring  nodes  in  the  finite  element  analysis.  The  parameters 
of  the  probe  are  set  to  mimic  Siemens  5-10MHz  probes.  The  probe  frequency  is 
7.27MHz,  the  sampling  rate  is  40MHz  and  the  fractional  bandwidth  is  60%. 

The  first  phantom  undergoes  uniform  compressions  in  the  z  direction  to 
achieve  strain  levels  of  2%  to  14%  in  2%  intervals.  Ground  truth  integer  dis¬ 
placement  values  are  used  as  the  initial  estimate  for  AM  to  decouple  the  perfor¬ 
mance  of  DP  from  AM.  Accurate  subpixel  displacement  field  is  calculated  with 
AM  and  the  mean  strain  values  are  compared  with  the  ground  truth  (Figure  1 
(a)-(c)).  The  results  are  only  shown  for  2%,  4%,  8%  and  14%  compression  for 
better  visualization.  The  results  with  two  threshold  values  for  IRLS  and  without 
IRLS  demonstrate  that  outlier  rejection  does  not  affect  the  mean  strain  value, 
while  increasing  the  regularization  weight  aa  increases  underestimation  of  the 
displacement.  The  rate  of  increase  of  the  underestimation  with  increasing  aa  is 
significantly  more  with  the  unbiased  regularization  (dashed  line)  as  expected. 

Significantly  higher  signal  to  noise  ratio  (SNR)  [2]  values  can  be  achieved 
with  outlier  rejection  (Figure  1  (d)-(f))  without  over-smoothing  the  image  with 
high  aa  values.  To  show  the  performance  of  the  overall  method,  the  initial  inte¬ 
ger  displacement  field  is  calculated  with  DP  and  accurate  displacement  field  is 
calculated  with  (Figure  1  (g)-(i)).  The  SNR  values  are  less  than  previous  case 
especially  at  high  strain  values,  where  DP  results  deviates  from  ground  truth. 

The  second  simulation  experiment  is  designed  to  show  the  effect  of  smooth¬ 
ness  weight  and  IRLS  threshold  on  contrast  to  noise  ratio  (CNR)  [2]  when  the 
correlation  is  lower  in  parts  of  the  image  due  to  fluid  motion.  The  phantom 
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(a)  T  =  0.005  (b)  T  =  0.01  (c)  no  IRLS 


(d)  T  =  0.005  (e)  T  =  0.01  (f)  no  IRLS 


(g)  T  =  0.005  (h)  T  =  0.01  (i)  no  IRLS 

Fig.  1.  Mean  and  SNR  of  the  elastograms  of  the  Field  II  simulated  uniform  phantom 
at  four  different  compression  levels  (shown  in  percentage)  for  three  IRLS  T  values.  The 
solid  and  dashed  lines  correspond  to  biased  and  unbiased  regularizations  respectively, 
(a) -(c)  shows  the  relative  underestimation  of  the  strain,  e  is  the  mean  strain  calculated 
with  the  elastography  method  and  e*  is  the  ground  truth,  (d)-(f)  shows  the  SNR  of 
the  AM.  (g)-(i)  shows  the  SNR  of  the  AM  with  initial  displacements  found  by  DP. 


contains  a  vein  oriented  perpendicular  to  the  image  plane  (Figure  2).  The  initial 
integer  displacement  is  generated  with  DP.  The  background  window  for  CNR 
calculation  is  located  close  to  the  target  window  to  show  how  fast  the  strain  is 
allowed  to  vary,  a  property  related  to  the  spatial  resolution.  The  maximum  CNR 
with  IRLS  is  5.3  generated  at  T  =  0.005  and  aa  =  38,  and  without  IRLS  is  4.8 
at  aa  =  338.  Such  high  aa  value  makes  the  share  of  the  data  term  in  the  cost 
function  very  small  and  causes  over-smoothing. 

Phantom  Results.  We  perform  freehand  palpation  experiment  on  a  breast 
phantom  to  examine  the  performance  of  the  frame  selection  technique.  50  frames 
of  RF  ultrasound  data  are  acquired  using  a  Siemens  Antares  system  (Issaquah, 
WA).  Our  custom  data  acquisition  program  is  connected  to  the  Axius  Direct 
Research  Interface  to  send  the  command  for  capturing  RF  data.  At  the  same 
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(a)  simulation  phan-  (b)  finite  element  strain 
tom 


(c)  CNR 


Fig.  2.  The  target  (circular)  and  background  (rectangular)  windows  for  CNR  calcula¬ 
tion  of  (c)  are  shown  in  (b) 


Fig.  3.  The  SNR  and  CNR  of  the  phantom  experiment  with  and  without  frame  selec¬ 
tion 


10  20 
width  (mm) 

(a)  CT  image  (b)  B-mode  image  (c)  strain  image 

Fig.  4.  Patient  experiment  results.  The  arrow  points  to  the  lumpectomy  cavity. 


time,  the  program  collects  tracking  information  from  a  Polaris  tracker  (Waterloo, 
Canada) .  Currently,  the  RF  frames  are  stored  on  the  ultrasound  system  and  are 
processed  offline.  Figure  3  shows  the  SNR  and  CNR  results.  In  automatic  frame 
selection,  Qij  (equation  8)  for  any  two  frames  in  a  buffer  of  size  15  frames 
is  calculated.  For  the  two  frames  which  give  the  minimum  Q,  the  strain  image 
is  obtained.  The  next  image  is  then  fed  to  the  buffer,  its  first  image  is  removed 
and  the  frame  selection  is  performed  again.  The  automatic  frame  selection  gives 
8  frame  pairs  for  strain  calculation  (as  seen  in  the  figure  by  8  SNR  and  CNR 
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values).  Without  frame  selection,  49  strain  images  are  calculated.  The  average 
CNR  and  SNR  values  are  improved  from  4.91  to  7.19  and  from  4.98  to  5.88  with 
frame  selection. 

Patient  Results.  We  have  acquired  freehand  palpation  ultrasound  RF  data 
using  the  Siemens  Antares  system  from  patients  approximately  four  weeks  after 
lumpectomy.  The  ultrasound  probe  is  tracked  with  the  Polaris  tracking  system. 
Optimal  frame  selection  is  performed  to  select  images  for  elastography  using  the 
AM  method.  The  strain  image  (Figure  4)  shows  that  the  AM  method  can  detect 
the  thin  hard  scar  tissue  even  though  it  is  close  to  the  cavity  fluids  which  undergo 
incoherent  motions  and  cause  signal  decorrelation.  Since  the  AM  method  finds 
the  displacement  of  all  the  samples  on  an  A-line  at  the  same  time,  the  correlated 
data  at  the  top  and  bottom  of  the  cavity  guide  the  method  to  find  the  correct 
displacement  inside  the  cavity  where  the  data  is  decorrelated. 

4  Discussion  and  Conclusion 

We  introduced  a  novel  method  for  calculating  a  dense  displacement  map  by 
analytic  minimization  of  a  cost  function.  We  used  the  IRLS  method  from  ro¬ 
bust  statistics  to  make  the  tracking  resistant  to  outliers.  Moreover,  we  exploited 
the  tracking  data  to  optimize  frame  selection.  Through  simulation  studies  using 
Field  II  and  finite  element  analysis,  we  showed  that  the  proposed  AM  method 
generates  high  quality  displacement  estimates.  The  elastography  method  works 
in  real-time.  A  comparison  of  the  IRLS  method  with  quality  guided  displacement 
tracking  [14]  which  also  aims  for  robustness  is  a  subject  of  future  work. 

We  chose  the  novel  application  of  the  lumpectomy  cavity  localization  as  the 
hard  scar  tissue  is  relatively  thin  and  demands  a  high  resolution  elastography 
method.  Also,  incoherent  fluid  motions  in  the  cavity  causes  large  decorrelations, 
requiring  a  robust  method.  We  have  an  active  Institutional  Review  Board(IRB) 
protocol  and  have  promising  results  from  9  patients  which  will  be  published  in 
future  work. 
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